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Abstract 

The  Coastal  Modeling  System  (CMS)  is  an  integrated  numerical  modeling 
system  for  simulating  nearshore  waves,  currents,  water  levels,  salinity  and 
sediment  transport,  and  morphology  change.  The  CMS  was  designed  and 
developed  for  coastal  inlets  and  navigation  applications,  including  channel 
performance  and  sediment  exchange  between  inlets  and  adjacent  beaches. 
The  present  report  provides  an  updated  description  of  the  mathematical 
formulations  and  numerical  methods  of  hydrodynamic,  salinity  and 
sediment  transport,  and  morphology  change  model  CMS-Flow.  The  CMS- 
Flow  uses  the  Finite  Volume  Method  on  Cartesian  grids  and  has  both  fully 
explicit  and  fully  implicit  time-stepping  schemes.  A  detailed  description  of 
the  explicit  time-stepping  scheme  was  provided  in  Militello  et  al.  (2004) 
and  Buttolph  et  al.  (2006).  The  present  report  focuses  on  the  recent 
changes  in  the  mathematical  formulations  and  the  implicit  time-stepping 
schemes.  The  CMS-Wave  and  CMS-Flow  models  are  tightly  coupled  within 
a  single  inline  code.  The  CMS-Wave  and  CMS-Flow  grids  may  be  the  same 
or  have  different  spatial  extents  and  resolutions.  The  hydrodynamic  model 
includes  physical  processes  such  as  advection,  turbulent  mixing,  combined 
wave-current  bottom  friction;  wave  mass  flux;  wind,  atmospheric 
pressure,  wave,  river,  and  tidal  forcing;  Coriolis  force;  and  the  influence  of 
coastal  structures.  The  implicit  hydrodynamic  model  is  coupled  to  a 
nonequilibrium  transport  model  of  multiple-sized  total-load  sediments. 
The  model  includes  physical  processes  such  as  hiding  and  exposure,  bed 
sorting  and  gradation,  bed  slope  effects,  nonerodible  surfaces,  and 
avalanching. 
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Preface 

This  report  was  written  by  the  Coastal  Inlets  Research  Program  (CIRP) 
and  funded  by  the  US  Army  Corps  of  Engineers,  Headquarters 
(HQUSACE).  The  CIRP  is  administered  at  the  US  Army  Engineer  Research 
and  Development  Center  (ERDC),  Coastal  and  Hydraulics  Laboratory 
(CHL)  under  the  Navigation  R&D  Program.  Jeff  E.  Mckee  was  HQUSACE 
Navigation  Business  Line  Manager  overseeing  CIRP.  Jeff  Lillycrop,  CHL, 
was  the  Technical  Director  of  the  Navigation  R&D  Program.  Dr.  Julie 
Rosati,  CHL,  was  the  CIRP  Program  Manager. 

CIRP  conducts  applied  research  to  improve  USACE  capabilities  to  manage 
federally-maintained  inlets  and  navigation  channels,  which  are  present  on 
all  coasts  of  the  United  States,  including  the  Atlantic  Ocean,  Gulf  of 
Mexico,  Pacific  Ocean,  Great  Lakes,  and  US  territories.  The  objectives  of 
CIRP  are  to  advance  knowledge  and  provide  quantitative,  predictive  tools 
to  accomplish  the  following:  (a)  support  management  of  federal  coastal 
inlet  navigation  projects  —  principally  the  design,  maintenance,  and 
operation  of  channels  and  jetties,  and  reduce  the  cost  of  dredging;  and  (b) 
preserve  the  adjacent  beaches  and  estuary  in  a  systems  approach  that 
treats  the  inlet,  beaches,  and  estuary  as  sediment-sharing  components.  To 
achieve  these  objectives,  CIRP  is  organized  in  work  units  conducting 
research  and  development  in  computational  modeling,  laboratory  and 
field  investigations,  and  technology  transfer. 

For  the  mission-specific  requirements,  CIRP  has  developed  a  finite- 
volume  model  based  on  nonlinear  shallow-water  flow  equations,  called 
CMS-Flow,  specifically  for  inlets,  navigation,  and  nearshore  project 
applications.  The  governing  equations  are  solved  using  both  explicit  and 
implicit  Finite  Volume  Method  on  Cartesian  coordinates.  The  model  is 
part  of  the  Coastal  Modeling  System  (CMS)  intended  to  simulate 
nearshore  waves,  flow,  sediment  transport,  and  morphology  change 
affecting  planning,  design,  maintenance,  and  reliability  of  federal 
navigation  projects. 

The  work  described  in  the  report  was  performed  under  the  Harbors 
Entrances  and  Structures  Branch  of  the  Navigation  Division,  US  Army 
Engineer  Research  and  Development  Center,  Coastal  and  Hydraulics 
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Laboratory  (ERDC-CHL).  At  the  time  of  publication,  Dr.  Donald  Ward  was 
Acting  Branch  Chief;  Dr.  Jackie  Pettway  was  Division  Chief;  Jeff  Lillycrop 
was  Technical  Director;  Dr.  Richard  Styles  was  Deputy  Director;  and  Jose 
Sanchez  was  the  Laboratory  Director. 

COL  Jeffrey  R.  Eckstein  was  ERDC  Commander.  Dr.  Jeffery  Holland  was 
ERDC  Director. 
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1  Introduction 

The  Coastal  Modeling  System  (CMS)  is  a  numerical  modeling  system  for 
nearshore  waves,  currents,  water  levels,  sediment  transport,  and 
morphology  change  (Militello  et  al.  2004;  Buttolph  et  al.  2006;  Lin  et  al. 
2008;  Reed  et  al.  2011).  The  system  was  developed  and  continues  to  be 
supported  by  the  Coastal  Inlets  Research  Program  (CIRP),  a  research  and 
development  program  of  the  US  Army  Corps  of  Engineers  (USACE)  that  is 
funded  by  the  Operation  and  Maintenance  Navigation  Business  Line  of  the 
USACE.  The  CMS  is  designed  for  coastal  inlets  and  navigation  applications 
including  channel  performance  and  sediment  exchange  between  inlets  and 
adjacent  beaches.  Modeling  provides  planners  and  engineers  with  essential 
information  for  reducing  costs  of  USACE  Operation  and  Maintenance 
activities.  CIRP  is  developing,  testing,  improving,  and  transferring  the  CMS 
to  Corps  Districts  and  industry  and  assisting  users  with  engineering  studies. 

The  overall  framework  of  the  CMS  and  its  components  are  presented  in 
Figure  1-1.  The  CMS  includes  a  flow  model  (CMS-Flow)  that  calculates 
hydrodynamics,  salinity  and  sediment  transport,  and  morphology  change 
and  a  spectral  wave  transformation  model  (CMS-Wave).  The  CMS-Wave 
and  CMS-Flow  models  are  tightly  coupled  within  a  single  inline  code  and 
may  have  the  same  or  different  computational  grids.  CMS  takes  advantage 
of  the  Surface-water  Modeling  System  (SMS)  interface  versions  8.1  and 
newer  (Militello  et  al.  2004;  Zundel  2006).  The  SMS  is  used  for  grid 
generation,  model  setup,  generating  input  files,  plotting  and  post¬ 
processing  of  results.  The  SMS  may  be  used  to  extract  boundary  condi¬ 
tions  from  a  larger-domain  simulation  (Buttolph  et  al.  2006).  The  SMS 
also  provides  a  link  between  the  CMS  and  the  Lagrangian  Particle 
Tracking  Model  (PTM)  (MacDonald  et  al.  2006;  Li  et  al.  2011). 

CMS-Wave  is  a  spectral  wave  transformation  model  and  solves  the  wave- 
action  balance  equation  using  a  forward  marching  Finite  Difference 
Method  (Mase  et  al.  2005;  Lin  et  al.  2008;  Lin  et  al.  2011a;  Lin  et  al. 

2012).  CMS-Wave  includes  physical  processes  such  as  wave  shoaling, 
refraction,  diffraction,  reflection,  wave-current  interaction,  wave  breaking, 
wind  wave  generation,  white  capping  of  waves,  and  the  influence  of  coastal 
structures  on  waves. 
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Figure  1-1.  CMS  framework  and  its  components. 


SMS 


CMS-Flow  is  a  two-dimensional  horizontal  (2DH)  depth-integrated  and 
wave-averaged  nearshore  hydrodynamic,  salinity  and  sediment  transport, 
and  morphology  change  model.  CMS-Flow  calculates  currents  and  water 
levels,  including  physical  processes  such  as  advection,  turbulent  mixing, 
combined  wave-current  bottom  friction;  wave  mass  flux;  wind,  atmospheric 
pressure,  wave,  river,  and  tidal  forcing;  Coriolis  force;  and  the  influence  of 
coastal  structures  (Buttolph  et  al.  2006;  Wu  et  al.  2011). 

CMS-Flow  has  three  noncohesive  sediment  transport  models  which  differ 
mainly  in  the  assumption  of  the  local  equilibrium  transport  for  the  bed  and 
suspended  loads  (Buttolph  et  al.  2006;  Sanchez  and  Wu  2011a, b).  CMS- 
Flow  can  simulate  any  number  of  sediment  size  fractions,  the  interactions 
between  size  fractions,  bed  sorting  and  layering,  and  morphology  change. 
The  sediment  transport  model  also  includes  processes  such  as  avalanching, 
nonerodible  surfaces,  and  bed  slope  effects. 

Typical  applications  of  CMS-Flow  include  analyses  of  past  and  future 
navigation  channel  performance;  wave,  current,  and  wave-current 
interaction  in  channels  and  in  the  vicinity  of  navigation  structures;  and 
sediment  management  issues  around  coastal  inlets  and  adjacent  beaches. 
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Some  examples  of  CMS-Flow  applications  are  as  follows:  Batten  and  Kraus 
(2006),  Buttolph  et  al.  (2006),  Zarillo  and  Brehin  (2007),  Li  et  al.  (2009), 
Beck  and  Kraus  (2010),  Byrnes  et  al.  (2010),  Dabees  and  Moore  (2011), 
Reed  and  Lin  (2011),  Rosati  et  al.  (2011),  Wang  and  Beck  (2011),  Beck  and 
Legault  (2012),  and  Lin  et  al.  (2013). 

A  detailed  description  of  the  theory  and  numerical  methods  for  CMS  was 
provided  in  Militello  et  al.  (2004)  and  Buttolph  et  al.  (2006).  Lin  et  al. 
(2008)  provided  a  description  of  the  mathematical  formulation  and 
verification  and  validation  cases  for  the  spectral  wave  model  CMS-Wave. 
More  recently,  several  verification  and  validation  tests  were  presented  for 
waves  (Lin  et  al.  2011b),  hydrodynamics  (Sanchez  et  al.  2011a),  and 
sediment  transport  and  morphology  change  (Sanchez  et  al.  2011b).  User 
guidance  on  how  to  setup,  run,  calibrate,  and  analyze  results  for  the  CMS 
was  presented  in  Lin  et  al.  (2011b)  and  Sanchez  et  al.  (2011a, b).  Since  the 
publication  of  Buttolph  et  al.  (2006),  there  have  been  several  changes  and 
added  features  and  processes  in  the  hydrodynamic  and  sediment  transport 
models.  The  purpose  of  the  present  report  is  to  provide  an  updated 
description  of  the  mathematical  formulation  and  numerical  methods  with 
focus  on  the  aspects  which  have  recently  changed.  In  particular,  this  report 
focuses  on  the  implicit  temporal  solution  scheme  for  hydrodynamics, 
sediment,  and  salinity  transport.  The  explicit  temporal  solution  scheme  is 
described  in  detail  in  Buttolph  et  al.  (2006).  Currently,  the  explicit  and 
implicit  versions  of  the  CMS-Flow  are  being  merged  into  a  single  code.  Once 
this  task  is  completed,  another  verification  and  validation  study  will  be 
conducted  with  comparisons  between  the  two  temporal  schemes  using  the 
previous  and  new  test  cases. 

To  put  both  the  purpose  and  contents  of  this  report  in  perspective,  it  is 
necessary  to  emphasize  that  this  report  is  not  a  self-standing,  all-inclusive 
document  about  the  implicit  CMS-Flow  model.  This  report  must  be 
employed  together  with  the  User’s  Guide  and  reports  published  for  the 
explicit  flow  model.  The  User’s  Guide  is  in  preparation  and  will  provide 
users  detailed  information  about  model  input  parameters  and  coefficients 
which  appear  in  this  report.  The  User’s  Guide  also  provides  guidance  on 
choice  of  methods,  formulas,  scaling  factors,  etc.  Last,  it  is  noted  that  CMS 
solves  a  depth-averaged  salinity  transport  equation  but  does  not  include 
an  equation  of  state  to  calculate  the  spatially  and  temporally  variable  water 
density.  These  cautions  are  provided  here  to  help  users  understand  model 
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limitations,  sources  of  information  needed  for  model  input,  and  proper 
application  of  the  model. 

The  report  is  organized  as  follows:  Chapter  2  presents  the  governing 
equations,  parameterizations,  empirical  equations,  boundary  and  initial 
conditions  for  hydrodynamics,  salinity  and  sediment  transport,  and 
morphology  change.  Variable  definitions  with  units  are  also  provided. 
Chapter  3  presents  the  numerical  methods.  The  Finite  Volume 
discretization  is  presented  for  a  general  transport  equation  in  order  to 
avoid  redundancy.  The  fully  implicit  iterative  solvers  for  hydrodynamics 
and  total-load  sediment  transport  are  described  in  detail.  Chapter  4  is  the 
summary,  and  Chapter  5  contains  the  references. 
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2  Mathematical  Formulations 

Coordinate  System 

Before  presenting  the  depth-integrated  and  wave-averaged  governing 
equations,  it  is  useful  to  define  the  coordinate  system  and  basic  variables. 
Variables  are  defined  spatially  in  a  Cartesian  coordinate  system 
x{  =  x  =  ( x,y,z ) ,  where  x  and  y  are  the  horizontal  coordinates,  and  z  is 

the  vertical  coordinate  (positive  is  upwards).  A  schematic  of  several  of  the 
main  variables  in  the  vertical  direction  is  provided  in  Figure  2-1.  The 
vertical  coordinate  datum  is  usually  the  Still  Water  Level  (SWL).  The  bed 
elevation  ( zb )  is  measured  from  the  vertical  datum  (i.e.,  negative 

downwards). 


Figure  2-1.  Vertical  conventions  used  for  the  bed  and  mean 
water  surface  elevation. 


Hydrodynamics 

This  section  provides  a  detailed  description  of  the  depth-integrated  and 
wave-averaged  hydrodynamic  equations.  The  equations  described  here 
closely  follow  those  presented  in  Phillips  (1977),  Mei  (1983),  and  Svendsen 
(2006).  Variable  definitions  and  the  final  hydrodynamic  equations  are 
provided  here.  Detailed  derivations  may  be  obtained  from  the  preceding 
references. 

Variable  Definitions 

The  instantaneous  current  velocity  ( ui )  is  split  into 
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U,  =  U- 


it; 


U; 


(2-1) 


in  which 

lij  =  current  (wave-averaged)  velocity  [m/s] 

«,  =  wave  (oscillatory)  velocity  [m/s]  with  wave-average  ut  =  o 
below  the  wave  trough 

u'  =  turbulent  velocity  fluctuation  [m/s]  with  ensemble  average 
(il  }  =  o,  and  wave  average  u'  =  o. 

The  wave-averaged  total  volume  flux  is  defined  as 


where: 

h  =  wave-averaged  total  water  depth  h  =  ff-zb  (Figure  2-1)  [m] 

Vi  =  total  mean  mass  flux  velocity,  or  simply  total  flux  velocity 
[m/s] 

rj  =  instantaneous  water  level  with  respect  to  the  Still  Water  Level 
(SWL)  [m] 

rj  =  wave-averaged  water  surface  elevation  with  respect  to  the  SWL 
(Figure  2-1)  [m] 

zh  =  bed  elevation  with  respect  to  the  SWL  (Figure  2-1)  [m]. 

The  total  flux  velocity  is  also  referred  to  as  the  mean  transport  velocity 
(Phillips  1977)  and  mass  transport  velocity  (Mei  1983).  The  current 
volume  flux  is  defined  as 


(2-3) 


where  Ut  is  the  depth-averaged  current  velocity.  Similarly,  the  wave 
volume  flux  is  defined  as 
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where: 

Uwi  =  depth-averaged  wave  flux  velocity  [m/s] 

77,  =  instantaneous  wave  trough  elevation  [m]  . 

Therefore,  the  total  flux  velocity  can  be  written  as 

Vl=Ui+Uwi.  (2-5) 


Governing  Equations 

On  the  basis  of  the  above  definitions,  and  assuming  depth-uniform 
currents,  the  general  depth-integrated  and  wave-averaged  continuity  and 
momentum  equations  can  be  written  as  (Phillips  1977;  Svendsen  2006) 


dh  d(hVj)_sm 
dt  dxj 


(2-6) 


dihVj)  d(  hVjVJ 

dt  dXj 


h  dpa 
P  (  )x, 


d 


dxj 


vth 


dv 

1  d 

Ox . 

J  / 

pdxj 

Sg  +  Rn  —p^Up^  +  —  -mb  — 
'  p  p 


(2-7) 


where: 

1  =  time  [s] 

Xj  =  horizontal  Cartesian  coordinate  in  th ejth  direction  [m],y=i,  2 

or  x,  y 

Sm  =  source  term  due  to  precipitation,  evaporation  and  structures 
(e.g.,  culverts)  [m/s] 

fc  =  2Q  sin  =  Coriolis  parameter  [rad/s]  in  which  Q  =  7.29x10-5 
rad/s  is  the  Earth’s  angular  velocity  of  rotation,  and  cj>  is  the 
latitude  in  degrees 

g  =  gravitational  constant  (-9.81  m/s2) 
pa  =  atmospheric  pressure  [Pa] 
p  =  water  density  (-1025  kg/ms) 
vt  =  total  eddy  viscosity  [m2/s] 
rsi  =  wind  surface  stress  [Pa] 


ERDC/CHL  TR-14-2 


8 


Sy  =  wave  radiation  stress  [Pa] 

Ry  =  surface  roller  stress  [Pa] 
mb  =  bed  slope  coefficient  [-] 

vhj  =  combined  wave  and  current  mean  bed  shear  stress  [Pa]. 

The  above  2DH  equations  are  similar  to  those  derived  by  Svendsen 
(2006),  except  for  the  inclusion  of  the  water  source/sink  term  in  the 
continuity  equation  and  the  atmospheric  pressure  and  surface  roller  terms 
and  the  bed  slope  coefficient  in  the  momentum  equation.  It’s  also  noted 
that  the  horizontal  mixing  term  is  formulated  slightly  differently  as  a 
function  of  the  total  flux  velocity,  similar  to  the  Generalized  Lagrangian 
Mean  (GLM)  approach  (Andrews  and  McIntyre  1978;  Walstra  et  al.  2000). 
This  approach  is  arguably  more  physically  meaningful  and  also  simplifies 
the  discretization  in  the  case  where  the  total  flux  velocity  is  used  as  the 
model  prognostic  variable. 

Bed  Shear  Stresses 

Bed  Roughness 

The  bed  roughness  is  specified  for  the  hydrodynamic  calculations  with 
either  a  Manning's  roughness  coefficient  ( n ),  Nikuradse  roughness  height 
(  ks ),  or  bed  friction  coefficient  ( ch ).  It  is  important  to  note  that  the  bed 

roughness  is  assumed  constant  in  time  and  not  changed  according  to  bed 
composition  and  bedforms.  This  is  a  common  engineering  approach  which 
can  be  justified  by  the  lack  of  data  to  initialize  the  bed  composition  and  the 
large  error  in  estimating  the  bed  composition  evolution  and  bedforms.  In 
addition  using  a  constant  bottom  roughness  simplifies  the  model  calibra¬ 
tion.  In  future  versions  of  CMS,  the  option  to  automatically  estimate  the  bed 
roughness  from  the  bed  composition  and  bedforms  will  be  added.  In 
addition,  the  bed  roughness  used  for  hydrodynamics  may  not  be  the  same 
as  that  which  is  used  for  the  sediment  transport  calculations  because  each 
sediment  transport  formula  was  developed  and  calibrated  using  specific 
methods  for  estimating  bed  shear  stresses  or  velocities,  and  these  cannot  be 
easily  changed. 

The  bed  friction  coefficient  ( cb )  is  related  to  the  Manning’s  roughness 
coefficient  ( n )  by  (Soulsby  1997) 

cb  =  gn2lrli3 


(2-8) 
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Commonly,  the  bed  friction  coefficient  is  calculated  by  assuming  a 
logarithmic  velocity  profile  as  (Graf  and  Altinakar  1998) 

I2 

c,  =  - - -  (2-9) 

6  [ln(z0  /h)  +  l\ 

where  k  =0.4  is  Von  Karman  constant,  and  z0  is  the  bed  roughness  length 
which  is  related  to  the  Nikuradse  roughness  ( ks )  by  z0  =  ks  /  30 
(hydraulically  rough  flow). 

Current-Related  Shear  Stress 

The  current  bed  shear  stress  is  given  by 

T(:l=pcbUU,  (2-10) 


where: 

p  =  water  density  (-1025  kg/m3) 
ch  =  bed  friction  coefficient  [-] 

U  =  \jU~U~  =  current  magnitude  [m/s]. 


The  magnitude  of  the  current-related  bed  shear  stress  is  simply 


G  =  Pcbu2 


(2-11) 


Wave-Related  Bed  Shear  Stress 

The  wave-related  bed  shear  stress  amplitude  is  given  by  (Jonsson  1966) 

Tw  =  \pfuK  (2-12) 

where  fw  is  the  wave  friction  factor,  and  uw  is  an  equivalent  or  representa¬ 
tive  bottom  wave  orbital  velocity  amplitude.  The  wave  friction  factor  ( fw )  is 
estimated  using  one  of  the  following: 


fw  =  exp(5.5r  °'2 


6.3 j  (Nielsen  1992) 


(2-13) 
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fw  =  0.237r~°'52  (Soulsby  1997) 


(2-14) 


fw  = 


exp(5.21r 
0.3 


-0.19 


6.0  forr>1.57 


forr  <  1.57 


(Swart  1974)  (2-15) 


where: 

r  =  relative  roughness  =  Aw  /  ks  [-] 
ks  =  Nikuradse  roughness  [m] 

Aw  =  semi-orbital  excursion  =  uwT/(2tt)  [m] 

T  =  wave  period  [s]. 

Mean  Bed  Shear  Stress  Due  to  Waves  and  Currents 

Under  combined  waves  and  currents,  the  mean  (wave-averaged)  bed  shear 
stress  is  enhanced  compared  to  the  case  of  currents  only.  This  enhancement 
of  the  bed  shear  stress  is  due  to  the  nonlinear  interaction  between  waves 
and  currents  in  the  bottom  boundary  layer.  In  CMS,  the  mean  (short-wave 
averaged)  bed  shear  stress  ( rbi )  is  calculated  as 

Tbi  =  Kc  Tci  (2-16) 


where: 

2  =  nonlinear  bottom  friction  enhancement  factor  ( 2  >  1 )  [-] 

wc  v  we  s  L-  -i 

rci  =  current-related  bed  shear  stress  [Pa]. 

The  nonlinear  bottom  friction  enhancement  factor  ( Awc )  is  calculated 

using  one  of  the  following  formulations  (name  abbreviations  are  given  in 
parenthesis): 

1.  Wu  et  al.  (20io)quadratic  formula  (QUAD) 

2.  Soulsby  (1995)  empirical  two-coefficient  data  fit  (DATA2) 

3.  Soulsby  (1995)  empirical  thirteen-coefficient  data  fit  (DATA13) 

4.  Fredsoe  (1984)  analytical  wave-current  boundary  layer  model  (F84) 

5.  Huynh-Thanh  and  Temperville  (1991)  numerical  wave-current  boundary 
layer  model  (HT91) 
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6. 

7- 

In 


Davies  et  al.  (1988)  numerical  wave-current  boundary  layer  model 
(DSK88) 

Grant  and  Madsen  (1979)  analytical  wave-current  boundary  layer  model 

(GM79) 

the  case  of  the  QUAD  formula,  Awc  is  given  by 


A, 


U 


(2-17) 


where  cw  is  an  empirical  coefficient,  and  uw  is  the  bottom  wave  orbital 
velocity  amplitude  based  on  linear  wave  theory.  For  random  waves, 
uw  ~  uws  where  uws  is  the  bottom  wave  orbital  velocity  amplitude  calculated 
based  on  the  significant  wave  height  and  peak  wave  period  (Equation  2-22). 
Wu  et  al.  (2010)  originally  proposed  setting  cw  =  0.5.  Here,  the  coefficient 

cw  has  been  calibrated  equal  to  1.33  for  regular  waves  and  0.65  for  random 
waves  to  agree  better  with  DATA2  formula. 


A  formula  similar  to  Equation  (2-17)  was  independently  proposed  by 
Wright  and  Thompson  (1983)  and  calibrated  using  field  measurements  by 
Feddersen  et  al.  (2000).  The  main  difference  in  the  two  formulations  is 
that  Wu  et  al.  (2010)  uses  the  bottom  wave  orbital  velocity  based  on  the 
significant  wave  height,  while  the  Wright  and  Thompson  (1983) 
formulation  uses  the  standard  deviation  of  the  bottom  orbital  velocity. 

The  DATA2,  DATA13,  F84,  HT91,  DSK88,  and  GM79  formulations  are 
calculated  using  the  general  parameterization  of  Soulsby  (1993): 

Awc  =  l  +  bXp(l-X)q  (2-18) 

where  X  =  tc/(tc  +  tw)  ,  and  b  ,  P ,  and  q  are  coefficients  given  by  (Soulsby 
et  al.  1993) 


X 


X±  +  X2  |cos <p\j)  +  (x3+  X4  |cos  ( p\r ) i°g10 


f 

J  w 


(2-19) 


where  x  =  ( b,p,q )  =  /(XiX'i’XsXa)  are  coefficients  which  have  been  fitted  to 
each  model  (Table  2-1). 


ERDC/CHL  TR-14-2 


12 


Table  2-1.  Fitting  coefficients  for  combined  wave-current  mean  bottom 

friction. 


Coefficient 

DATA2 

DATA13 

F84 

HT91 

DSK88 

GM79 

bi 

1.2 

0.47 

0.29 

0.27 

0.22 

0.73 

b2 

0.0 

0.69 

0.55 

0.51 

0.73 

0.40 

bi 

0.0 

-0.09 

-0.10 

-0.10 

-0.05 

-0.23 

b4 

0.0 

-0.08 

-0.14 

-0.24 

-0.35 

-0.24 

Pi 

0.0 

-0.53 

-0.77 

-0.75 

-0.86 

-0.68 

P2 

0.0 

0.47 

0.10 

0.13 

0.26 

0.13 

P3 

0.0 

0.07 

0.27 

0.12 

0.34 

0.24 

P4 

0.0 

-0.02 

0.14 

0.02 

-0.07 

-0.07 

<7/ 

3.2 

2.34 

0.91 

0.89 

-0.89 

1.04 

q: 

0.0 

-2.41 

0.25 

0.40 

2.33 

-0.56 

qs 

0.0 

0.45 

0.50 

0.50 

2.60 

0.34 

qi 

0.0 

-0.61 

0.45 

-0.28 

-2.50 

-0.27 

j 

0.0 

8.8 

3.00 

2.70 

2.70 

0.50 

The  GM79,  DATA2,  and  DATA13  models  use  the  logarithmic  relationship 
for  the  bed  friction  coefficient  given  by  Equation  (2-9).  In  the  case  of  the 
F84,  HT91,  and  DSK88  models,  the  bed  friction  coefficient  is  linearly 
interpolated  in  log-space  using  the  tabulated  values  presented  in  Soulsby 
(1997). 


In  the  case  of  the  F84,  HT91,  DSK88,  and  GM79  models,  the  wave  friction 
factors  are  linearly  interpolated  in  log-space  using  the  tabulated  values 
found  in  Soulsby  (1997).  In  the  case  of  the  DATA2  and  DATA13  formulas, 
the  wave  friction  factor  is  estimated  using  Equation  (2-13). 


Bottom  Wave  Orbital  Velocity 

The  bottom  wave  orbital  velocity  amplitude  for  regular  waves  (uw)  is 
calculated  based  on  linear  wave  theory  as 


jiH 

T  sinh  {kh) 


(2-20) 


where: 


H  =  wave  height  [m] 
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T  =  wave  period  [s] 
k  =  wave  number  [rad/m]. 

Unless  specified  otherwise,  for  random  waves,  uw  is  set  to  an  equivalent  or 
representative  bottom  orbital  velocity  amplitude  equal  to  uw  =  V2 unm 
where  urms  the  root-mean-squared  bottom  wave  orbital  velocity  amplitude 
is  defined  here  following  Soulsby  (1987;  1997): 

_  r*oo 

urms  =  var(ub)  =  £  Su(f)df  (2-21) 


where: 

var()  =  variance  function, 

uh  =  instantaneous  bottom  orbital  velocity  [m/s] 

Su  =  wave  orbital  velocity  spectrum  density  [s  m2/s2] 

/  =  wave  frequency  [l/s] . 

It  is  noted  that  the  definition  of  urms  is  slightly  different  from  others  such  as 

Madsen  (1994),  Myrhaug  et  al.  (2001),  and  Wiberg  and  Sherwood  (2008) 
which  include  factor  of  2  in  their  definition.  A  simple  approximation  for  urms 

from  linear  wave  theory  and  the  root-mean-squared  wave  height 
Hrms  =  Hs  /  V2  (for  a  Rayleigh  distribution)  is  given  by 


u 


JiH,. 


T  V2  sinh(/ch) 


(2-22) 


Wiberg  and  Sherwood  (2008)  reported  that  urms  estimates  using  Hrms  and 
Tp  agree  reasonably  well  with  field  measurements  (except  for  Tp  <8.8  s) 
and  produces  better  estimates  than  other  combinations  with  Hrms ,  Hs ,  Tp , 
and  the  zero-crossing  wave  period  (71).  The  zero-crossing  wave  period  is 

calculated  as  the  average  period  (time  lapse)  between  consecutive  upward 
or  downward  intersections  of  the  water  level  time  series  with  the  zero  water 
line.  A  better  approach  is  to  assume  a  spectral  shape  such  as  the  Joint  North 
Sea  Wave  Project  (JONSWAP)  (Hasselman  et  al.  1973)  and  obtain  an 
explicit  curve  for  urms  by  summing  the  contributions  from  each  frequency 
(Soulsby  1987;  Wiberg  and  Sherwood  2008).  A  simple  explicit  expression  is 
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provided  below  based  on  the  JONSWAP  (y  -  3.3)  spectrum  following  the 
work  of  Soulsby  (1987): 


n  H, 

Urms  =0.134-^ 


1  +  tanh 


-7.76— ^-  +  1. 34 


(2-23) 


where  Tn  =  ^|hJg  .  The  above  expression  agrees  closely  with  the  curves 
presented  by  Soulsby  (1987;  1997). 

In  some  cases  the  bottom  wave  orbital  velocity  amplitude  is  calculated 
based  on  the  significant  wave  height  and  peak  wave  period  ( uws )  as 


*HS 

Tp  sinh  (kh) 


(2-24) 


Bed  slope  Friction  Coefficient 

It  is  noted  that  in  the  presence  of  a  sloping  bed,  the  bottom  friction  acts  on 
a  larger  surface  area  for  the  same  horizontal  area.  This  increase  in  bottom 
friction  is  included  through  the  coefficient  (Mei  1989;  Wu  2007) 


mb  =  |  Vzb 


+ 


dz. 


dy 


+  1 


(2-25) 


where  zb  is  the  bed  elevation,  and  V 


'  d_  d_ 
dx ’ dy ’ 


For  bottom  slopes  of 


1/5  and  1/3,  the  above  expression  leads  to  an  increase  in  bottom  friction  of 
2.0  percent  and  5.4  percent,  respectively.  In  most  morphodynamic 
models,  the  bottom  slope  is  assumed  to  be  small,  and  the  above  term  is 
neglected.  However,  it  is  included  here  for  completeness. 


Eddy  Viscosity 

The  term  eddy  viscosity  arises  from  the  fact  that  small-scale  vortices  or 
eddies  on  the  order  of  the  grid  cell  size  are  not  resolved,  and  only  the 
large-scale  flow  is  simulated.  The  eddy  viscosity  is  intended  to  simulate  the 
dissipation  of  energy  at  smaller  scales  than  the  model  can  simulate.  In  the 
nearshore  environment,  large  mixing  or  turbulence  occurs  due  to  waves, 
wind,  bottom  shear,  and  strong  horizontal  gradients.  Therefore,  the  eddy 
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viscosity  is  an  important  parameter  which  can  have  a  large  influence  on 
the  calculated  flow  field  and  resulting  sediment  transport.  In  CMS-Flow, 
the  total  eddy  viscosity  ( vt )  is  equal  to  the  sum  of  three  parts:  l)  a  base 

value  ( v0 );  2)  the  current-related  eddy  viscosity  ( vc );  and  3)  the  wave- 
related  eddy  viscosity  ( vw )  defined  as  follows: 


vt=yo+vc+vw  (2-26) 

The  base  value  ( v0 )  is  approximately  equal  to  the  kinematic  viscosity 

(~i.8ixicr6  m2/s)  but  may  be  changed  by  the  user.  The  other  two 
components  (vc  and  vw)  are  described  in  the  sections  below. 


Current-Related  Eddy  Viscosity  Component 

There  are  four  algebraic  models  for  the  current-related  eddy  viscosity:  1) 
Falconer  Equation;  2)  depth-averaged  parabolic;  3)  subgrid;  and  4)  mixing- 
length.  The  default  turbulence  model  is  the  subgrid  model  but  may  be 
changed  by  the  user. 

Falconer  Equation 

The  Falconer  (1980)  equation  was  the  default  method  applied  in  earlier 
versions  of  CMS  (Militello  et  al.  2004)  for  the  current-related  eddy 
viscosity.  The  equation  is  given  by 

vc  =  0.575 cbUh  (0-07 


where  cb  is  the  bottom  friction  coefficient,  U  is  the  depth-averaged  current 
velocity  magnitude,  and  h  is  the  total  water  depth. 


Depth-averaged  Parabolic  Model 

The  second  option  for  the  current-related  eddy  viscosity  is  the  depth- 
averaged  parabolic  model  given  by 


vc  =  cvmch 


(2-28) 


where  u,c  =  ^ tc  /  p  is  the  bed  shear  velocity,  and  cv  is  approximately 

equal  to  k/6=o.o667  but  is  set  as  a  calibrated  parameter  whose  value  can 
be  up  to  1.0  in  irregular  waterways  with  weak  meanders  or  even  larger  for 
strongly  curved  waterways. 
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Subgrid  Model 

The  third  option  for  calculating  the  current-related  eddy  viscosity  ( vc )  is 
the  subgrid  turbulence  model  given  by 


v=c,mh+(chAf 


(2-29) 


where: 


cv  =  vertical  shear  coefficient  [-] 


c'<  =  horizontal  shear  coefficient  [-] 

A  =  (average)  grid  size  [m] 

1%^ 

e  =  deformation  (strain  rate)  tensor  = 


lidV.  dV 
- L  + - L 

2  dXj  dxi 


The  empirical  coefficients  cv  and  ch  are  related  to  the  turbulence  produced 
by  the  bed  shear  and  horizontal  velocity  gradients.  The  parameter  cv  is 
approximately  equal  to  k/6=o.o667  (default)  but  may  vary  from  o.oi  to 
0.2.  The  variable  ch  is  equal  to  approximately  the  Smagorinsky  coefficient 
(Smagorinsky  1963)  and  may  vary  between  0.1  and  0.3  (default  is  0.2). 

Mixing  Length  Model 

The  mixing  length  model  implemented  in  CMS  for  the  current-related 
eddy  viscosity  includes  a  component  due  to  the  vertical  shear  and  is  given 
by  (Wu  2007) 


(2-30) 


where: 

^  =  mixing  length  ( =  k  min (chh,yr) )  [m] 
y'  =  distance  to  the  nearest  wall  [m] 
ch  =  horizontal  shear  coefficient  [-]. 

The  empirical  coefficient  ch  is  usually  between  0.3  and  1.2.  The  effects  of 
bed  shear  and  horizontal  velocity  gradients,  respectively,  are  taken  into 
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account  through  the  first  and  second  terms  on  the  right-hand  side  of 
Equation  (2-30).  It  has  been  found  that  the  modified  mixing  length  model 
is  better  than  the  depth-averaged  parabolic  eddy  viscosity  model  that 
accounts  for  only  the  bed  shear  effect. 

Wave-Related  Eddy  Viscosity 

The  wave  component  of  the  eddy  viscosity  is  separated  into  two 
components: 


vw=cwfuwsHs+cbrh 


D 


br 


1/3 


(2-31) 


where: 

cwf  =  wave  bottom  friction  coefficient  for  eddy  viscosity  [-] 
uws  =  peak  bottom  orbital  velocity  [m/s]  based  on  the  significant 
wave  height  Hs  [m]  and  peak  wave  period  T  [s] 
cbr  =  wave  breaking  coefficient  for  eddy  viscosity  [-] 

Dbr  =  wave  breaking  dissipation  [N/m/s]. 


The  first  term  on  the  right-hand  side  of  Equation  (2-31)  represents  the 
component  due  to  wave  bottom  friction,  and  the  second  term  represents  the 
component  due  to  wave  breaking.  The  coefficient  cwf  is  approximately  equal 

to  0.5  and  may  vary  from  0.5  to  2.0.  The  coefficient  cbr  is  approximately 
equal  to  0.1  and  may  vary  from  0.04  to  0.15. 


Wave  Radiation  Stresses 

The  wave  radiation  stresses  ( St] )  are  calculated  using  linear  wave  theory  as 
(Longuet-Higgins  and  Stewart  1961;  Dean  and  Dalrymple  1984) 


s„  =  /iXt/.o 


1 

U9  2 


dfde 


(2-32) 


where: 

/  =  the  wave  frequency  [l/s] 
0  =  the  wave  direction  [rad] 
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K 

Hs 

w. 


n 


g 


c 

k 


wave  energy  =  ^pgH2  [N/m] 
significant  wave  height  [m] 
wave  unit  vector  =  (cos  6,  sin  9)  [-] 
J 1  for  i  =  j 
[o  for  i  j 

=  2kh  |  [-] 
c  2{  sinh2  kh) 

wave  group  velocity  [m/s] 

wave  celerity  [m/s] 
wave  number  [rad/m]. 


The  wave  radiation  stresses  and  their  gradients  are  computed  within  the 
wave  model  and  interpolated  in  space  and  time  in  the  flow  model. 

Roller  Stresses 

As  a  wave  transitions  from  non-breaking  to  fully-breaking,  some  of  the 
energy  is  converted  into  momentum  that  goes  into  the  aerated  region  of 
the  water  column.  This  phenomenon  is  known  as  the  surface  roller.  The 
surface  roller  contribution  to  the  wave  stresses  ( Rtj )  is  given  by 


Rij=2Esr™i™j 


(2-33) 


where: 

Esr  =  surface  roller  energy  density  [N/m] 

Wj  =  wave  unit  vector  =  (cos  6m ,  sin  9m )  [-]. 

The  surface  roller  is  calculated  within  CMS-Wave.  An  effect  of  the  surface 
roller  is  to  shift  the  peak  alongshore  current  velocity  closer  to  shore. 
Another  side  effect  of  the  surface  roller  is  to  improve  model  stability 
(Sanchez  et  al.  2011a).  The  influence  of  the  surface  roller  on  the  mean 
water  surface  elevation  is  relatively  minor  (Sanchez  et  al.  2011a). 

Wave  Flux  Velocity 

In  the  presence  of  waves,  the  oscillatory  wave  motion  produces  a  net  time- 
averaged  mass  (volume)  transport  referred  to  as  Stokes  drift.  In  the  surf 
zone,  the  surface  roller  also  provides  a  contribution  to  the  mean  wave 
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mass  flux.  The  mean  wave  mass  flux  velocity,  or  simply  the  mass  flux 
velocity,  is  defined  as  the  mean  wave  volume  flux  divided  by  the  local 
water  depth  and  is  approximated  here  as  (Phillips  1977;  Svendsen  2006) 


£7  ,  = 


(Ew+2Esr)wi 

phc 


(2-34) 


where: 

Ew  =  wave  energy  =^pgH;  [N/m] 

Hs  =  significant  wave  height  [m] 

Esr  =  surface  roller  energy  density  [N/m] 
wt  =  (cos  #,sin  6)  =  wave  unit  vector  [-] 
c  =  wave  celerity  [m/s]. 

The  first  component  is  due  to  the  Stokes  velocity  while  the  second 
component  is  due  to  the  surface  roller  (only  present  in  the  surf  zone). 

Wind  Surface  Stress 

The  wind  surface  stress  is  calculated  as 

Tsi  =  PaCDWW{  (2-35) 


where: 

pa  =  air  density  at  sea  level  [-1.2  kg/m3] 

CD  =  wind  drag  coefficient  [-] 

Wf  =  10-m  wind  speed  [m/s] 

W  =  \jw~W~  =  10-m  wind  velocity  magnitude  [m/s]. 

The  wind  speed  is  calculated  using  either  an  Eulerian  or  Lagrangian 
reference  frame  as 


w;  =wiE  -Ywui 


(2-36) 


where: 
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WfE  =  io-m  atmospheric  wind  speed  relative  to  the  solid  earth 
(Eulerian  wind  speed)  [m/s] 

yw  =  equal  to  zero  for  the  Eulerian  reference  frame  or  one  for  the 

Lagrangian  reference  frame 
Ui  =  current  velocity  [m/s]. 


Using  the  Lagrangian  reference  frame  or  relative  wind  speed  is  more 
accurate  and  realistic  for  field  applications  (Bye  1985;  Pacanowski  1987; 
Dawe  and  Thompson  2006),  but  the  option  to  use  the  Eulerian  wind  speed 
is  provided  for  idealized  cases  and  backward  compatibility.  The  drag 
coefficient  is  calculated  using  the  formula  of  Hsu  (1988)  and  modified  for 
high  wind  speeds  based  on  field  data  by  Powell  et  al.  (2003)  (Figure  2-2): 


C 


D 


(  \2 

K 


(14.56  — 2ln  PE  J 


for  PE  <  30  m/s 


10  3max(3.86-0.04PE,i.5)  forPE>30m/s 


(2-37) 


Figure  2-2.  Modified  Hsu  (1988)  wind  drag  coefficient. 


Powell  et  al.  (2003)  speculate  that  the  reason  for  the  decrease  in  drag 
coefficient  with  higher  wind  speeds  is  due  to  increasing  foam  coverage 
leading  to  the  formation  of  a  slip  surface  at  the  air-sea  interface. 

Wind  measurements  taken  at  heights  other  than  10  m  are  converted  to  to¬ 
rn  wind  speeds  using  the  1/7  rule  (Shore  Protection  Manual  (SPM)  1984; 
Coastal  Engineering  Manual  (CEM)  2002): 


l  z  J 


PE  =  W; 


(2-38) 
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where  z  is  the  elevation  above  the  sea  surface  of  the  wind  measurement, 
and  Wr  is  the  wind  velocity  at  height  z  . 

Boundary  Conditions 

Wall  Boundary 

The  wall  boundary  condition  is  a  closed  boundary  and  is  applied  at  any  cell 
face  between  wet  and  dry  cells.  Any  unassigned  boundary  cell  at  the  edge  of 
the  model  domain  is  assumed  to  be  closed  and  is  assigned  a  wall  boundary. 
A  zero  normal  flux  to  the  boundary  is  applied  at  closed  boundaries.  Two 
boundary  conditions  are  available  for  the  tangential  flow: 

1.  Free-slip:  no  tangential  shear  stress  (wall  friction) 

2.  Partial-slip:  tangential  shear  stress  (wall  friction)  calculated  based  on  the 
log  law. 

Assuming  a  log  law  for  a  rough  wall,  the  partial-slip  tangential  shear  stress 
is  given  by 


Twan=f>cwauui  (2-39) 

where  t/,  is  the  magnitude  of  the  wall  parallel  current  velocity,  and  cwall  is 
the  wall  friction  coefficient  equal  to 


Cwall  ~ 


K 


In  (yP/y0) 


(2-40) 


Here,  v0  is  the  roughness  length  of  the  wall  and  is  assumed  to  be  equal  to 
that  of  the  bed  (i.e.,  v0  =  z0 ).  The  distance  from  the  wall  to  the  cell  center  is 

yp- 


Flux  Boundary 

The  flux  boundary  condition  is  typically  applied  to  the  upstream  end  of  a 
river  or  stream  and  specified  as  either  a  constant  or  time  series  of  total 
water  volume  flux  ( Q ).  In  a  2DH  model,  the  total  volume  flux  needs  to  be 

distributed  across  the  boundary  in  order  to  estimate  the  depth-averaged 
velocities.  This  is  done  using  a  conveyance  approach  in  which  the  current 
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velocity  is  assumed  to  be  related  to  the  local  flow  depth  ( h )  and  Manning’s 
(n)  as  U  cch'  In .  Here,  r  is  an  empirical  conveyance  coefficient  equal  to 

approximately  2/3  for  uniform  flow.  The  smaller  the  r  value  the  more 
uniform  the  current  velocities  are  across  the  flux  boundary.  The  water 
volume  flux  (g, )  at  each  boundary  cell  (/')  is  calculated  as 


qt  =  hut 


•f RampQ 


h] 


r+ 1 


Ec 


hr+1 

en)— — AZ 
n. 


n,. 


(2-41) 


where: 

i  =  subscript  indicating  a  boundary  cell 
qt  =  volume  discharge  at  boundary  cell  i  per  unit  width  [m2/s] 

e  =  unit  vector  for  inflow  direction  =  (sin  rp,  cos  (p) 

(p  =  inflow  direction  measured  clockwise  from  North  [deg] 

n  =  boundary  face  unit  vector  (positive  outward) 

Q  =  specified  total  volume  flux  across  the  boundary  [m3/s] 

n  =  Manning’s  coefficient  [s/mhs] 
r  =  empirical  constant  equal  to  approximately  2/3 
A /  =  cell  width  in  the  transverse  direction  normal  to  flow  [m] 
fRamp  =  ramp  function  [-]  (described  in  Chapter  3). 

The  total  volume  flux  is  positive  into  the  computational  domain.  Since  it  is 
not  always  possible  to  orient  all  flux  boundaries  to  be  normal  to  the  inflow 
direction,  the  option  is  given  to  specify  an  inflow  direction  (<p).  The  angle 

is  specified  in  degrees  clockwise  from  true  North.  If  the  angle  is  not 
specified,  then  the  inflow  angle  is  assumed  to  be  normal  to  the  boundary. 
The  total  volume  flux  is  conserved  independently  of  the  inflow  direction. 

Water  Level  Boundary 

The  general  formula  for  the  boundary  water  surface  elevation  is  given  by 

Vb  =  fRamp^E  +  ^  ~  fRampM 0  (2-42) 


where: 


tjB  =  boundary  water  surface  elevation  [m] 
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rfE  =  specified  external  boundary  water  surface  elevation  [m] 

AT)  =  water  surface  elevation  offset  [m] 
ff0  =  initial  boundary  water  surface  elevation  [m] 
rjc  =  correction  to  the  boundary  water  surface  elevation  based  on 
the  wind  and  wave  forcing  [m] 

rjG  =  water  surface  elevation  component  derived  from  user  specified 
gradients  [m] 

f i{amP  -  ramp  function  [-]  (described  in  Chapter  3). 

The  external  water  surface  elevation  (ffE)  may  be  specified  as  a  time  series, 

both  spatially  constant  and  varying  or  calculated  from  tidal/harmonic 
constituents.  When  a  time  series  is  specified,  the  values  are  interpolated 
using  piecewise  Lagrangian  polynomials.  By  default,  second  order 
interpolation  is  used  but  can  be  changed  by  the  user.  If  tidal  constituents 
are  specified,  then  rfE  is  calculated  as 


(2-43) 


where: 


i  =  subscript  indicating  a  tidal  constituent 
4  =  mean  amplitude  [m] 

4  =  node  (nodal)  factor  [-] 
coi  =  frequency  [deg/hr] 

t  =  elapsed  time  from  midnight  of  the  starting  year  [hrs] 

Vt°  +  ui  =  equilibrium  phase  [deg] 

=  phase  lag  [deg]. 

The  mean  amplitude  and  phase  may  be  specified  by  the  user  or  interpolated 
from  a  tidal  constituent  database.  The  nodal  factor  is  a  time-varying  correc¬ 
tion  to  the  mean  amplitude.  The  equilibrium  phase  has  a  uniform  com¬ 
ponent  ( V° )  and  a  relatively  smaller  periodic  component.  The  zero- 

superscript  of  V°  indicates  that  the  constituent  phase  is  at  time  zero. 

Table  2-2  provides  a  list  of  tidal  constituents  currently  supported  in  CMS. 
More  information  on  US  tidal  constituent  values  can  be  obtained  from  the 
US  National  Oceanographic  and  Atmospheric  Administration  (NOAA) 
(http://tidesoniine.nos.noaa.gov)  and  National  Ocean  Service  (NOS)  (htto://co- 
oos.nos.noaa.gov). 
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Table  2-2.  Tidal  constituents  implemented  in  CMS.  Constituent  speeds  are  in  degrees  per  mean  solar  hour. 


Constituent 

Speed 

Constituent 

Speed 

Constituent 

Speed 

Constituent 

Speed 

SA* 

0.041067 

SSA* 

0.082137 

MM* 

0.54438 

MSF* 

1.0159 

MF* 

1.098 

2Q1* 

12.8543 

Ql* 

13.3987 

RH01* 

13.4715 

01* 

13.943 

Ml* 

14.4967 

PI* 

14.9589 

SI* 

15.0 

15.0411 

Jl* 

15.5854 

001* 

16.1391 

2N2* 

27.8954 

MU2* 

27.9682 

N2* 

28.4397 

NU2* 

28.5126 

M2 

28.9841 

LDA2* 

29.4556 

L2* 

29.5285 

* 

CM 

1— 

29.9589 

S2 

30 

R2* 

30.0411 

K2 

30.0821 

2SM2* 

31.0159 

2MK3* 

42.9271 

M3* 

43.4762 

MK3* 

44.0252 

MN4* 

57.4238 

M4 

57.9682 

MS4* 

58.9841 

84* 

60.0 

M6 

86.9523 

S6* 

90.0 

M8* 

115.9364 

If  a  harmonic  boundary  condition  is  applied,  then  the  node  factors  are  set 
to  one  and  the  equilibrium  arguments  are  set  to  zero.  The  harmonic 
boundary  condition  is  provided  as  an  option  for  simulating  idealized  or 
hypothetical  conditions. 

The  water  surface  elevation  offset  ( A  77 )  is  assumed  spatially  and  temporally 

constant  and  may  be  used  to  correct  the  boundary  water  surface  elevation 
for  vertical  datums  and  sea  level  rise.  The  component  rjG  is  intended  to 

represent  regional  gradients  in  the  water  surface  elevation,  is  assumed  to  be 
constant  in  time,  and  is  only  applicable  when  ffE  is  spatially  constant.  When 

applying  a  water  level  boundary  condition  to  the  nearshore,  local  flow 
reversals  and  boundary  problems  may  result  if  the  wave-  and  wind-induced 
setup  is  not  included.  This  problem  is  avoided  by  adding  a  correction  ( rfc ) 

to  the  local  water  level  to  account  for  the  wind  and  wave  setup  as 

pghB  +  VG  +  Vc )  =  ^  +  Twx  ~  T,,  (2-44) 

where  h„  is  the  boundary  total  water  depth,  and  rsx ,  twx  ,  and  rhx  are  the 
wind,  wave,  and  bottom  stresses  in  the  boundary  direction  (x).  The  wave 
forcing  term  is  equal  to 


n 


dxj 


Sii+R,-phUwiUWJ 


(2-45) 
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The  correction  rjc  is  only  applicable  when  tje  is  spatially  constant  as  in  the 
case  of  a  single  water  surface  elevation  time  series. 


Cross-shore  Boundary 

In  the  implicit  flow  solver,  a  cross-shore  boundary  condition  is  applied  in 
the  nearshore  by  solving  the  lD  cross-shore  momentum  equation 
including  wave  and  wind  forcing  (Wu  et  al.  2010).  Along  a  cross-shore 
boundary,  it  is  assumed  that  a  well-developed  longshore  current  exists 
(quasi-steady  conditions  with  longshore  gradients  in  advection,  diffusion, 
and  water  levels  equal  to  zero).  These  assumptions  are  valid  for  relatively 
long  coasts  with  shore-parallel  contours  and  simplify  the  alongshore  (y- 
direction)  momentum  equation  to 


d 

dx 


pv,hdV‘ 


dx 


Tsy  +  TWy  ~  Tby 


(2-46) 


where  rsy ,  rWy ,  and  xby  are  the  surface,  wave,  and  bottom  stresses  in  the 

longshore  direction,  respectively.  Equation  (2-46)  is  solved  iteratively  to 
determine  the  longshore  current  velocity.  The  cross-shore  (x)  component  of 
the  velocity  is  assigned  a  zero-gradient  boundary  condition.  The  longshore 
current  velocity  is  applied  when  the  flow  is  directed  inwards.  When  the  flow 
is  directed  outwards,  a  zero-gradient  boundary  condition  is  applied  to  the 
longshore  current  velocity. 


The  water  level  due  to  waves  and  wind  at  the  cross-shore  boundary  can  be 
determined  by  assuming  a  zero  alongshore  gradient  of  flow  velocity  and 
negligible  cross-shore  current  velocity.  For  this  case,  the  cross-shore 
momentum  equation  reduces  to 


7  drj 

P9hyr-  =  Ts 
dx 


'  TWx  Xbx 


(2-47) 


where  tsx  ,  zWx ,  and  zbx  are  the  surface,  wave,  and  bottom  stresses  in  the 
cross-shore  direction.  The  water  level  boundary  condition  is  applied  for 
the  case  when  the  flow  is  directed  outwards.  When  the  flow  is  directed 
inwards,  a  zero-gradient  boundary  condition  is  applied  to  the  water  level 
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Salinity  Transport 

Overview 

The  characteristics  of  salinity  are  important  in  the  coastal  environment 
because  salinity  can  impact  marine  plants  and  animals  and  influence  the 
dynamic  behavior  of  cohesive  sediments.  Because  modifications  of  coastal 
inlets,  such  as  channel  deepening  and  widening  and  rehabilitation  or 
extension  of  coastal  structures,  may  alter  the  salinity  distribution  within 
estuaries  or  bays,  it  is  often  useful  and  convenient  to  simulate  the  salinity 
within  the  scope  of  an  engineering  project  to  determine  if  a  more  detailed 
water  quality  modeling  study  is  necessary.  It  is  important  to  emphasize 
that  the  CMS  is  not  intended  to  be  used  as  a  water  quality  model.  The  CMS 
solves  the  depth-averaged  (2DH)  salinity  transport  equation  and  should  be 
used  only  for  cases  where  the  water  column  is  well  mixed.  If  there  is  flow 
stratification,  a  3D  model  should  be  utilized.  It  is  also  noted  that  the 
salinity  is  not  used  to  update  the  water  density  which  is  assumed  to  be 
constant.  Thus  any  horizontal  water  density  gradients  due  to  varying 
salinity  on  the  hydrodynamics  are  assumed  to  be  negligible. 


Transport  Equation 

CMS  calculates  the  salinity  transport  based  on  the  following  2DH  salinity 
conservation  equation  (Li  et  al.  2012): 


d(hclal) ,  dfhVf^-) 


a 


dt 


dXj 


dxj 


VSalh°C^L 

dx.i 


+s, 


sal 


(2-48) 


where: 

=  depth-averaged  salinity  [usually  in  ppt  or  psu] 

h  =  water  depth  [m] 

Vj  =  total  flux  velocity  [m/s] 

vsal  =  horizontal  mixing  coefficient  vsal  =vtl  <jsal  [m2/s] 
vt  =  total  eddy  viscosity  [m2/s] 

osal  =  Schmidt  number  for  salinity  (approximately  equal  to  1.0)  [-] 
Ssal  =  source/sink  term  [ppt  m/s]. 

The  above  equation  represents  the  horizontal  fluxes  of  salt  in  water  bodies 
and  is  balanced  by  exchanges  of  salt  via  diffusive  fluxes.  Major  processes 
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that  influence  the  salinity  are  as  follows:  seawater  exchange  at  ocean 
boundaries,  freshwater  inflows  from  rivers,  precipitation  and  evaporation 
at  the  water  surface,  and  groundwater  fluxes  (not  included  here). 

Initial  Condition 

The  initial  condition  for  salinity  transport  may  be  specified  as  a  constant,  a 
spatially  variable  dataset  usually  calculated  from  a  previous  simulation  or 
by  solving  a  2DH  Laplace  equation: 


V2csa/  =  0  (2-49) 

where  V2  is  the  Laplace  operator.  The  equation  is  solved  given  any 
number  of  user-specified  initial  salinity  values  at  locations  within  the 
model  domain,  using  the  initial  salinity  values  at  open  boundaries  and 
applying  a  zero-gradient  boundary  condition  at  all  closed  boundaries. 

Boundary  Conditions 

At  cell  faces  between  wet  and  dry  cells,  a  zero-flux  boundary  condition  is 
applied.  A  salinity  time  series  must  be  specified  at  all  open  boundaries  and 
is  applied  when  the  flow  is  directed  inward  of  the  modeling  domain.  If  the 
flow  is  directed  outward  of  the  modeling  domain,  then  a  zero-gradient 
boundary  condition  is  applied. 

Sediment  Transport 

Overview 

For  sand  transport  in  coastal  waters,  the  wash-load  (i.e.,  sediment 
transport  which  does  not  appreciably  contribute  to  the  bed  material)  can 
be  assumed  to  be  negligible,  and,  therefore,  the  total -load  transport  is 
equal  to  the  sum  of  the  bed-  and  suspended-load  transports.  There  are 
currently  three  sediment  transport  models  available  in  CMS: 

1.  equilibrium  total  load  (ET) 

2.  equilibrium  bed  load  plus  non-equilibrium  suspended  load  (EBNS) 

3 .  non-equilibrium  total-load  (NET) . 

The  main  differences  among  the  three  models  are  the  assumption  of  local 
equilibrium  transport  for  the  bed  and  suspended  loads.  The  ET  model 
assumes  that  both  the  bed  and  suspended  loads  are  in  local  equilibrium 
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and  solves  a  simple  mass  conservation  equation  known  as  the  Exner 
equation.  It  is  the  least  computationally  intensive  but  can  lead  to  model 
instabilities  due  to  the  assumption  of  local  equilibrium  transport.  The 
EBNS  model  solves  a  transport  equation  for  the  depth-averaged 
suspended-load  concentration,  while  the  bed  load  is  assumed  to  be  in 
equilibrium  and  is  included  in  a  bed  change  equation.  Using  a  non- 
equilibrium  formulation  for  suspended  load  is  more  realistic  and  provides 
more  accurate  results.  However,  it  still  assumes  that  the  bed  load  is  in 
local  equilibrium.  A  more  realistic  approach  is  to  use  non-equilibrium 
formulations  for  both  the  bed  and  suspended  loads  as  in  the  NET  model. 
The  NET  implemented  in  the  CMS  combines  the  bed-  and  suspended-load 
transport  equations  into  a  single  total-load  transport  equation  thus 
reducing  the  computational  costs.  The  ET  and  EBNS  models  are  only 
available  using  a  single  sediment  size  and  with  the  explicit  hydrodynamic 
temporal  scheme.  They  are  described  in  detail  in  Buttolph  et  al.  (2006) 
and  are  not  repeated  here.  The  third  model  is  available  with  both  the 
explicit  and  implicit  schemes.  The  support  for  multiple  grain  sizes  is 
currently  only  in  the  NET  model  with  the  implicit  time-stepping  scheme. 
The  multiple-sized  non-equilibrium  sediment  transport  model  is 
introduced  in  this  section. 

Non-equilibrium  Total-Load  Transport  Model 

Total-load  Transport  Equation 

The  single-sized  sediment  transport  model  described  in  Sanchez  and  Wu 
(2011a)  was  extended  to  multiple-sized  sediments  within  CMS  by  Sanchez 
and  Wu  (2011b).  In  this  model,  the  sediment  transport  is  separated  into 
current-  and  wave-related  transports.  The  transport  due  to  currents 
includes  the  stirring  effect  of  waves,  and  the  wave-related  transport 
includes  the  transport  due  to  asymmetric  oscillatory  wave  motion  as  well 
as  steady  contributions  by  Stokes  drift,  surface  roller,  and  undertow.  The 
current-related  bed  and  suspended  transports  are  combined  into  a  single 
total -load  transport  equation,  thus  reducing  the  computational  costs  and 
simplifying  the  bed  change  computation.  The  2DH  transport  equation  for 
the  current-related  total  load  is 


d 

f  hCtk ) 

,  d(hUftk) 

d 

v  h  ^rsk^tk  ^ 

dt 

,  Ptk  , 
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for  j  =  1,2;  k  =  l,2,...,N ,  where  N is  the  number  of  sediment  size  classes 
and 


Ctk  =  actual  depth-averaged  total-load  sediment  concentration 

[kg/m3]  for  size  class  k  defined  as  Ctk  =  qtk  /  (Uh)  in  which  qtk 
is  the  total-load  mass  transport 

Ctk,  =  equilibrium  depth-averaged  total-load  sediment  concentration 
[kg/m3]  for  size  class  k  and  described  in  the  equilibrium 
concentration  and  transport  rates  section 
(3tk  =  total -load  correction  factor  described  in  the  Total-Load 

Correction  Factor  section  [-] 

rsk  =  fraction  of  suspended  load  in  total  load  for  size  class  k 

described  in  fraction  of  suspended  sediments  section  [-] 
vs  =  horizontal  sediment  mixing  coefficient  described  in  the 

horizontal  sediment  mixing  coefficient  section  [m2/s] 
a ,  =  total-load  adaptation  coefficient  described  in  the  adaptation 

coefficient  section  [-] 
cosk  =  sediment  fall  velocity  [m/s]. 

In  the  above  equation,  the  first  term  represents  the  temporal  variation  of 
Ctk ;  the  second  term  represents  the  horizontal  advection;  the  third  term 

represents  the  horizontal  diffusion  and  dispersion  of  suspended  sediments; 
and  the  last  term  represents  the  erosion  and  deposition.  The  equation  may 
be  applied  to  single-sized  sediment  transport  by  using  a  single  sediment  size 
class  (i.e.,  N  =  l).  The  bed  composition,  however,  does  not  vary  when  using  a 
single  sediment  size  class.  The  units  of  sediment  concentration  used  here 
are  kg/m3  rather  than  dimensionless  volume  concentrations  in  order  to 
avoid  numerical  precision  errors  at  low  concentrations. 

In  the  above  equations,  it  is  assumed  that  the  wave  flux  velocity  is  not 
included  in  the  momentum  equations.  If  the  wave  flux  velocity  is  included, 
then  the  total  flux  velocity  ( V )  should  be  used  instead  of  the  depth- 
averaged  current  velocity  ( U ).  The  reason  for  this  is  because  without  a 
wave-induced  sediment  transport  to  counter  the  offshore  directed 
transport  due  to  the  undertow,  the  model  would  predict  excessive 
movement  of  sediment  offshore. 
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Fraction  of  Suspended  Sediment 

In  order  to  solve  the  system  of  equations  for  sediment  transport  implicitly, 
the  fraction  of  suspended  sediment  must  be  determined  explicitly.  This  is 
done  by  assuming 


QSk  Qsk* 


1  sk  —  — 


Qtk  Qtk* 


(2-51) 


where  qsk  and  qtk  are  the  actual  suspended-  and  total -load  transport  rates 
and  qsk,  and  qtk,  are  the  equilibrium  suspended-  and  total -load  transport 
rates. 


Adaptation  Coefficient 

The  total -load  adaptation  coefficient  ( a  t)  is  an  important  parameter  in  the 
sediment  transport  model.  There  are  many  variations  of  this  parameter  in 
literature  (Lin  1984;  Gallappatti  and  Vreugdenhil  1985;  and  Armanini  and 
di  Silvio  1986).  CMS  uses  a  total -load  adaptation  coefficient  ( a, )  that  is 
related  to  the  total -load  adaptation  length  ( Lt )  and  time  ( Tt )  by 


Uh 


ut : 


at(os 


(2-52) 


where: 

cos  =  sediment  fall  velocity  corresponding  to  the  transport  grain  size 

for  single-sized  sediment  transport  or  the  median  grain  size 
for  multiple-sized  sediment  transport  [m/s] 

U  =  depth-averaged  current  velocity  [m/s] 
h  =  water  depth  [m]. 

The  adaptation  length  or  time  is  a  characteristic  distance  or  time  for 
sediment  to  adjust  from  non-equilibrium  to  equilibrium  transport. 

Because  the  total  load  is  a  combination  of  the  bed  and  suspended  loads, 
the  associated  adaptation  length  may  be  calculated  as  Lt  =  rsLs  +  ( 1  -  r  )Lb  or 

Lt  =  max(Zs ,  Lb ) ,  where  Ls  and  Lb  are  the  suspended-  and  bed-load 
adaptation  lengths.  The  symbol  Ls  is  defined  as 
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LS=  —  =  UTS  (2-53) 

acos 

in  which  a  and  Ts  are  the  adaptation  coefficient  lengths  for  suspended 
load.  The  adaptation  coefficient  ( a )  can  be  calculated  either  empirically 
or  based  on  analytical  solutions  to  the  pure  vertical  convection-diffusion 
equation  of  suspended  sediment.  One  example  of  an  empirical  formula  is 
that  proposed  by  Lin  (1984): 


a  =  3.25  +  0.55ln 


co 


s 


KU * 


(2-54) 


where  u,  is  the  bed  shear  stress,  and  k  is  the  von  Karman  constant. 
Armanini  and  di  Silvio  (1986)  proposed  an  analytical  equation: 


1 

a 


1.5 


-1/6 


(2-55) 


where  S  is  the  thickness  of  the  bottom  layer  defined  by  8  =  33z0 ,  and  z0  is 

the  zero-velocity  distance  from  the  bed.  Gallappatti  (1983)  proposed  the 
following  equation  to  determine  the  suspended-load  adaptation  time: 


T  = 

s 


h 

— exp 

u. 


(1.57  - 20.12u. ) a),3  +(326.832u, 22047  -0.2)o).2 
+  (0.1385 In ur  — 6.406l)«,  +(o.5467u  -2.1963) 


(2-56) 


where  u*  is  the  current  related  bottom  shear  velocity,  ur  =  uJU ,  and 
co*  =  cos  /  u* . 

The  bed-load  adaptation  length  ( Lb )  is  generally  related  to  the  dimension  of 

bed  forms  such  as  sand  dunes.  Large  bed  forms  are  generally  proportional 
to  the  water  depth,  and,  therefore,  the  bed-load  adaptation  length  can  be 
estimated  as  Lh  =  abh  in  which  ab  is  an  empirical  coefficient  on  the  order  of 

5-10.  Although  limited  guidance  exists  on  methods  to  estimate  Lh ,  the 
determination  of  Lh  is  still  empirical  and  in  the  developmental  stage.  For  a 

detailed  discussion  of  the  adaptation  length,  the  reader  is  referred  to  Wu 
(2007).  In  general,  it  is  recommended  that  the  adaptation  length  be 
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calibrated  with  field  data  in  order  to  achieve  the  best  and  most  reliable 
results. 

Total-Load  Correction  Factor 

The  total-load  correction  factor  ( j3tk )  accounts  for  the  vertical  distribution 

of  the  suspended  sediment  concentration  and  velocity  profiles  as  well  as 
the  fact  that  bed  load  travels  at  a  slower  velocity  than  the  depth-averaged 
current  velocity  (Figure  2-3)  (Wu  2007).  By  definition,  j3tk  is  the  ratio  of 

the  depth-averaged  total -load  and  flow  velocities. 


Figure  2-3.  Schematic  of  sediment  and  current 
vertical  profiles. 


In  a  combined  bed-load  and  suspended-load  model,  the  correction  factor 
is  given  by 


1 

^  rjpsk+(l ~rsk)U/ubk  (2  57 ' 

where  ubk  is  the  bed-load  velocity,  and  J3sk  is  the  suspended-load  correction 

factor  and  is  defined  as  the  ratio  of  the  depth-averaged  suspended  sediment 
and  flow  velocities.  Since  most  sediment  is  transported  near  the  bed,  both 
the  total-  and  suspended-load  correction  factors  ( f3tk  and  j8sk )  are  usually 

less  than  1  and  typically  in  the  range  of  0.3  and  0.7,  respectively.  By 
assuming  logarithmic  current  velocity  and  exponential  suspended  sediment 
concentration  profiles,  an  explicit  expression  for  the  suspended-load 
correction  factor  ( j3sk )  may  be  obtained  as  (Sanchez  and  Wu  2011b) 


■f„..1‘c‘<fe  E,(4/D-E,(^)  +  ln(/l/Z)e-«J-ln(l/Z)e 

uf’  c,dz  e-*a[ln(l/Z)-l][l-e-‘(1-«] 

J  zb+a 
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where: 

fa=(Dskh/ek  [-] 

A  =  a/h  [-] 

Z  =  zjh  [-] 

sk  =  vertical  mixing  coefficient  [m2/s] 

a  =  reference  height  near  the  bed  for  the  suspended  load  [m] 
zfl  =  apparent  roughness  length  [m] 

r*  oo  g  ^ 

Ei(x)=  I  — dt  (exponential  integral). 

J  X  f 

The  equation  can  be  further  simplified  by  assuming  that  the  reference 
height  is  proportional  to  the  apparent  roughness  length  (e.g.,  a  =  30 zQ )  so 

that  psk  =  psk  (Z,$ 4)  •  Figure  2-4  shows  a  comparison  of  the  suspended-load 

correction  factor  based  on  the  logarithmic  velocity  with  exponential  and 
Rouse  suspended  sediment  concentration  profiles.  Both  cases  show 
similar  behavior  for  the  suspended-load  correction  factor  ( j3sk ).  For  fine 

sediments  (small  fall  velocity),  J3sk  is  close  to  1.0  and  experiences  smaller 
influences  by  the  bottom  roughness,  while  for  course  sediments,  J3sk  can 
be  as  low  as  0.5  and  is  largely  influenced  by  the  bottom  roughness.  This  is 
to  be  expected  since  course  sediments  are  transported  more  closely  to  the 
bottom,  compared  to  fine  sediments. 


Figure  2-4.  Suspended  load  correction  factors  based  on  the  logarithmic  velocity  profile  and  (a) 
exponential  and  (b)  Rouse  suspended  sediment  profile.  The  Rouse  number  is  r  =  cos  /  {tcut) . 


The  bed-load  velocity  ( uhk )  in  Equation  (2-57)  is  calculated  using  the  van 
Rijn  (1984a)  formula  with  re-calibrated  coefficients  from  Wu  et  al.  (2006): 
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where: 


(2-59) 


5  =  specific  gravity  [-] 
g  =  gravitational  constant  (-9.81  m/s2) 

dk  =  characteristic  grain  diameter  for  the  kth  size  class  [m] 
zb  =  {n'  /  n)3/2  rb  =  grain-related  bed  shear  stress  [Pa] 
n' =  ^so6  /  20  =  grain-related  Manning’s  roughness  coefficient 
[s/mbs] 

zcrk  =  critical  bed  shear  stress  for  the  k,h  size  class  [Pa]. 


Bed  Change  Equation 

The  fractional  bed  change  is  calculated  as 


r)  c)y 

=  at ask  (C*  —  C*. )  +  Dsqhk  -^-k- 


(2-60) 


where: 

zh  =  bed  elevation  with  respect  to  the  vertical  datum  [m] 
p'm  =  bed  porosity  [-] 

ps  =  sediment  (material)  density  (-2650  kg/m3  for  quartz 
sediment) 

Ds  =  empirical  bed-slope  coefficient  (constant)  [-] 
qbk  =  hUCtk{±-rsk)  is  the  bed-load  mass  transport  rate  magnitude 
[kg/m/s]. 

The  first  term  on  the  right-hand  side  of  the  above  equation  represents  the 
bed  change  due  to  sediment  exchange  near  the  bed.  The  last  term  accounts 
for  the  effect  of  the  bed  slope  on  bed-load  transport.  The  bed-slope 
coefficient  ( D  )  is  usually  about  0.1  to  3.0.  For  a  detailed  derivation  of  the 

above  equation,  the  reader  is  referred  to  Sanchez  and  Wu  (2011a).  The  total 
bed  change  is  calculated  as  the  sum  of  the  bed  changes  for  all  size  classes: 


(2-61) 
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Bed  Material  Sorting  and  Layering 

Bed  sorting  is  the  process  in  which  the  bed  material  changes  size 
composition  (fraction  of  each  grain  size  class).  The  bed  is  discretized  into 
multiple  layers  to  consider  the  heterogeneity  of  bed  material  size 
composition  along  the  bed  depth.  The  fraction  of  each  size  class  is  then 
calculated  and  stored  in  each  layer.  The  sorting  of  sediments  is  then 
calculated  using  the  mixing  or  active  layer  concept  (Hirano  1971;  Karim 
and  Kennedy  1982;  Wu  1991).  The  mixing  layer  is  the  top  layer  of  the  bed 
which  exchanges  directly  with  the  sediment  moving  in  the  water  column. 
In  other  words,  only  the  sediment  in  the  mixing  layer  exchanges  with  the 
moving  sediment  in  the  water  column;  whereas,  the  sediment  in  the 
subsurface  layers  below  the  mixing  layer  does  not  directly  exchange  or 
contact  with  the  moving  sediment. 

The  temporal  variation  of  the  bed-material  gradation  in  the  first  (mixing 
or  active)  layer  is  calculated  as  (Wu  2007) 


where: 

St  =  thickness  of  the  first  layer  [m] 
plk  =  fraction  of  the  kth  sediment  size  in  the  first  layer  [-] 

\plk  for  dzb /dt-d61/dt>0 

Pik  =  [-] 

p2k  for  dzb  /  dt  —  d81  /  dt  <  0 

p2k  =  fraction  of  the  kth  sediment  size  in  the  second  layer  [-]. 

The  first  term  in  Equation  (2-62)  corresponds  to  the  erosion  and  deposition 
due  to  the  k,h  sediment  size  while  the  second  term  corresponds  to  the 
sediment  exchange  between  the  first  and  second  bed  layers  (Figure  2-5). 

At  the  beginning  of  each  time-step,  the  mixing  or  active  layer  thickness  is 
calculated  as 


(2-62) 


d(5i Pit)  _  ( dzb)  *  ggi  dzb 

dt  (dt)k  Plk  dt  dt 


81  =  min  [max(5min  ,2d50 ,0.5//,.  ),§max 


(2-63) 
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Figure  2-5.  Multiple  bed  layer  model  of  bed  material  sorting  (after  Wu  2007). 


where  Hr  is  the  ripple  height,  and  Smm  and  Smax  are  user-specified 
minimum  and  maximum  layer  thicknesses,  respectively. 


The  bed-material  sorting  in  the  second  layer  is  calculated  as 


d(S2p2k)  

* 

dS1  dzb 

dt 

Pik 

dt  dt 

(2-64) 


where  S2  is  the  thickness  of  the  second  layer.  It  is  noted  that  there  is  no 
material  exchanged  between  the  sediment  layers  below  the  second  layer. 

The  sediment  transport,  bed  change,  and  bed  gradation  equations  are 
solved  simultaneously  (coupled)  but  are  decoupled  from  the  flow 
calculation  at  each  time-step. 

Sediment  Fall  Velocity 

The  sediment  fall  velocity  may  be  user-specified  or  is  calculated  using  the 
formula  by  Soulsby  (1997): 


v 

“•  =  d 


1/2 

10. 362  +  1.049d,3  )  -10.36 


(2-65) 


where: 

v  =  kinematic  viscosity  [m2/s] 

d  =  grain  size  [m] 

d„  =  dimensionless  grain  size  [-]. 
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The  dimensionless  grain  size  is  given  by 


d, 


=  d 


(s-l)g 


1/3 


(2-1) 


where: 

^  =  sediment  specific  gravity  or  relative  density  [-] 
g  =  gravitational  constant  (9.81  m/s2). 

Equilibrium  Concentrations  and  Transport  Rates 

In  order  to  close  the  system  of  equations  describing  the  sediment 
transport,  bed  change,  and  bed  sorting  equations,  the  fractional 
equilibrium  depth-averaged  total -load  concentration  ( Ctk* )  must  be 

estimated  from  an  empirical  formula.  The  depth-averaged  equilibrium 
concentration  is  defined  as 


C,f  =  ^  (2-67) 

where  qtk„  is  the  total -load  transport  for  the  kth  sediment  size  class 
estimated  from  an  empirical  formula.  For  convenience,  Ctk,  is  written  in 
general  form  as 


Ctk*  —  Pik^tk  (2-68) 

where  p]k  is  the  fraction  of  the  sediment  size  (k)  in  the  first  (top)  bed  layer, 
and  C*.  is  the  potential  equilibrium  total -load  concentration.  The  potential 
concentration  ( C*k )  can  be  interpreted  as  the  equilibrium  concentration 
for  uniform  sediment  of  size  dk .  The  above  equation  is  essential  for  the 
coupling  of  sediment  transport,  bed  change,  and  bed  sorting  equations. 

Lund-CIRP 

Camenen  and  Larson  (2005,  2007,  2008)  developed  general  sediment 
transport  formulas  for  bed  and  suspended  loads  under  combined  waves 
and  currents.  These  are  referred  to  as  the  Lund-CIRP  transport  formulas. 
The  general  transport  formulas  can  be  used  for  both  symmetric  and 
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asymmetric  waves,  but  for  simplicity  the  waves  are  assumed  to  be 
symmetric.  The  current-related  bed-  and  suspended-load  transport  with 
wave  stirring  is  given  by 


chr 


^(s-l)gcl 


3 

50 


fbp,  12^0  cwm  exp 


0 

-4.5-  cr 


0 


CW  ) 


(2-69) 


q '  -  fsPfifl  — 

60,. 


^(s-l)gd 


3 

50 


1-exp 


co  h 


(2-70) 


where: 

qb*  =  equilibrium  bed-load  transport  rate  [kg/m/s] 

=  equilibrium  suspended-load  transport  rate  [kg/m/s] 

=  Shields  parameters  due  to  currents  [-] 

_  mean  shields  parameters  due  to  waves  and  currents  [-] 

0CW  =  maximum  Shields  parameters  due  to  waves  and  currents  [-] 

0cr  =  critical  Shields  parameter  [-] 

s  =  vertical  sediment  diffusivity  [m2/s] 
cR  =  reference  bed  concentration  [kg/m3] 

fb  =  bed-load  scaling  factor  (default  l.o)  [-] 

fs  =  suspended-load  scaling  factor  (default  l.o)  [-]. 


The  critical  Shields  parameter  is  calculated  using  Equation  (2-100).  The 
mean  and  maximum  Shields  parameters  are  calculated  as 

®cw,m  =  V0c+0lm+20c0u,,mcos  q>  (2-71) 

0ciy  =  V0c+0^+20c0u,cos  cp  (2-72) 

The  mean  wave  Shields  parameter  is  calculated  as  0U,  m  =0^/2  assuming 

a  sinusoidal  wave.  The  Shields  parameters  for  currents  and  waves  are 
given  by 


0 


■  c\w 


c\w 


g(ps-p)d 


(2-73) 
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in  which  the  subscript  c|  iv  indicates  either  the  current-  (c)  or  wave-related 
(in)  component.  The  current-related  shear  stress  (  tc  )  is  calculated  with 

Equation  (2-11).  The  wave-related  bed  shear  stress  is  calculated  with 
Equation  (2-12)  and  the  wave  friction  factor  ( fw )  of  Swart  (1974)  given  by 

Equation  (2-15). 

The  total  bed  roughness  is  assumed  to  be  a  linear  summation  of  the  grain- 
related  roughness  (ksg),  form-drag  (ripple)  roughness  (£„.),  and  sediment- 

related  roughness  ( kss ): 


K,c\ w  —  Kg  +  ^sr,c\w  +  Ks,c\w  (2"74) 

Here,  the  grain-related  roughness  is  estimated  as  ksg  =  2 d50 .  The  ripple 
roughness  ( ksr )  is  calculated  as  (Soulsby  1997) 


H2 

k  .  =  7.5^^- 


“V,c|u; 


v sr,c\w 


(2-75) 


where  Hr  and  Lr  are  the  ripple  height  and  length,  respectively. 

The  current-  and  wave-related  sediment  roughnesses  are  estimated  as 

Ks,c\W  =  5d50®  c\w  (2-76) 


The  above  equation  must  be  solved  simultaneously  with  the  expressions 
for  the  bottom  shear  stress  because  the  roughness  depends  on  the  stress. 

The  reference  concentration  is  given  by 


C*  =  Ar®c w,m  exp 


-4.5 


0„ 


(2-77) 


where  the  coefficient  AcR  is  determined  by  the  following  relationship: 

AcR  =0.0035  exp  (—0.3d.)  (2-78) 


The  vertical  sediment  diffusivity  is  calculated  as 
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£  =  h\De 


1/3 


where  De  is  the  total  effective  dissipation  given  by 


De  ~  k bDbr  +kcDc+  k^Dw 


(2-79) 


(2-80) 


in  which  kb  ,  kc ,  and  kw  are  coefficients;  Dbr  is  the  wave  breaking 
dissipation  (from  the  wave  model);  and  Dc  and  Dw  are  the  bottom  friction 

dissipation  due  to  currents  and  waves,  respectively.  The  dissipation  from 
bottom  friction  due  to  current  ( D  )  and  the  dissipation  from  bottom 

friction  due  to  waves  ( Dw )  are  expressed  as 


Dc\w  =  *c\wu*c\w  (2-81) 

The  coefficient  kb  =0.017  (Camenen  and  Larson  2008),  and  kc  and  kw  are 
a  function  of  the  Schmidt  number: 


k  —  —a 

c\w  _  Q  o\w 


(2-82) 


where  cr  is  either  the  current-  or  wave-related  Schmidt  number 


f  1 

c\w 


calculated  from  the  following  relationships  (Camenen  and  Larson  2008): 


°c|  w  — 


ac\w+bc\w  sin2 


71  CO„ 


2  u. 


c\w 


for  <  1 


u 


c\w 


(2-83) 


71  u 


c\w 


2 


for  — >  1 


u 


*c\w 


with  the  coefficients  ac  =  0.4 ,  bc  =  3.5 ,  aw  —  0.15  ,  and  bw  =  1.5 . 

For  multiple-sized  (non-uniform)  sediments,  the  fractional  equilibrium 
sediment  transport  rates  are  calculated  as  (Wu  and  Lin  2011) 


Qbk* 


Mk  'PikPs12^®^  exp 


0 

-4.5-  crk 


0 


(2-84) 
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Qsk* 


^(s-ljgd, 


fs&PlkPsCRkU- 


(0. 


sk 


1-exp 


™skh 


(2-85) 


where: 

k  =  subscript  indicating  the  sediment  size  class 
4  =  hiding  and  exposure  function  [-] 

rsk  =  fraction  of  suspended  load  for  each  size  class  defined  by 
Equation  (2-51)  [-] 

Plk  =  fraction  of  the  kth  sediment  size  in  the  first  layer  [-]. 

The  availability  of  sediment  fractions  is  included  through  p]k ,  while  hiding 
and  exposure  of  grain  sizes  are  accounted  for  by  directly  multiplying  the 
transport  rates. 

Van  Rijn 

The  van  Rijn  (1984a, b)  equations  for  bed-load  and  suspended-load 
transport  are  used  with  the  recalibrated  coefficients  of  van  Rijn  (2007a, b) 
as  given  by 


qb*  =  fbPs°-015Uh 


U  -U 

e  cr 

1.5 

Ko) 

yj(s-l)gd50 

h  J 

1.2 


(2-86) 


qs,  =  fn,0.012Ud: 


50 


£7  -£/ 


y](s  —  l)gd50 


2.4 


-0.6 


(2-87) 


where: 

Ucr  =  critical  depth-averaged  velocity  for  incipient  motion  [m/s] 

Ue  =  effective  depth  averaged  velocity  [m/s]. 

The  effective  depth-averaged  velocity  is  calculated  as  Ue=U  +  yuw  with 
7=0.4  for  random  waves  and  y  =  0.8  for  regular  waves.  The  bottom  wave 
orbital  velocity  based  on  linear  wave  theory  is  uw-  For  random  waves, 
uw  =  uws  where  uws  is  based  on  the  significant  wave  height  and  peak  wave 
period  (see  Equation  2-24).  The  critical  depth-averaged  velocity  is 
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estimated  as  Ucr  =  j3cUcrc  +  (1  -  J3c)ucrw  where  =U  /  (U  +  uw)  is  a 

blending  factor.  The  critical  depth-averaged  current  velocity  ( Ucrc )  is 
given  by  Equation  (2-102),  and  the  critical  bottom-wave-orbital  velocity 
amplitude  ( ucrw )  is  given  by  Equation  (2-103). 

According  to  van  Rijn  (2007a),  the  bed-load  transport  formula  predicts 
transport  rates  by  a  factor  of  2  for  velocities  higher  than  0.6  m/s,  but 
under  predicts  transport  rates  by  a  factor  of  2  to  3  for  velocities  close  to  the 
initiation  of  motion. 

The  van  Rijn  formulas  (1984a, b;  2007a, b)  were  originally  proposed  for  well- 
sorted  sediments.  The  sediment  availability  is  included  by  multiplication  of 
transport  rates  with  the  fraction  of  the  sediment  size  class  in  the  upper  bed 
layer.  The  hiding  and  exposure  are  considered  by  a  correction  factor  which 
multiples  to  the  critical  velocity.  When  applied  to  multiple-sized  sediments, 
the  fractional  equilibrium  transport  rates  are  calculated  as 


Qbk*  =  fbPsPik°-015Uh 


1.5 

ftU 

V( s-l)gdk 

UJ 

1.2 


(2-88) 


qsk*  =  fsPsPlk0.012Uh 


u,-&2ucrk 

2.4 

\dk) 

^(s-l)gdk 

[h) 

7-0.6 


(2-89) 


The  availability  of  sediment  fractions  is  included  through  p]k ,  while  hiding 

and  exposure  of  grain  sizes  are  accounted  for  by  multiplying  the  critical 
velocity  ( \Jcrk )  by  a  correction  function  ( ). 

Soulsby-van  Rijn 

Soulsby  (1997)  proposed  the  following  equation  for  the  total  load  sediment 
transport  rate  under  action  of  combined  current  and  waves: 


qb*  =  fbPs°-005Uh 


(  \ 

U  -U 

e  crc 

2.4 

\d50) 

-\/(s  i)gd50 

h  J 

qs*=fsPs0.012Uh 


U  -U 

e  crc 

2.4 

\d50) 

^s-l)gd50j 

h 

(2-90) 


(2-91) 
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where: 


Ue  =  jU2  +  °'°18u2rms  =  effective  velocity  [m/s] 

V  cb 

u,ms  =  root-mean-square  bottom  wave  orbital  [m/s] 

Ucrc  =  critical  depth-averaged  velocity  for  initiation  of  motion  for 
currents  based  on  Van  Rijn  (1984c)  [m/s]. 

The  bed  friction  coefficient  ( cb )  is  calculated  using  Equation  (2-9)  with  the 
bed  roughness  length  ( z0 )  set  to  0.006  m  following  Soulsby  (1997). 

The  Soulsby- van  Rijn  formula  is  modified  for  multiple-sized  sediments 
similarly  to  the  van  Rijn  formula  in  the  previous  section  with  the  equation 


qbk*=fbPsPik  0-005UH 


2.4 

\dk) 

yjs-l)gdk 

[h) 

1.2 


(2-92) 


<lsk*=fsPsPik0-012Uh 


\u,-eknucrk\ 

2.4 

\dk) 

js-l)gdk 

d-0.6 

u*k 


(2-93) 


As  in  the  case  of  the  van  Rijn  transport  formula,  the  availability  of 
sediment  fractions  is  included  through  p]k ,  while  hiding  and  exposure  of 

grain  sizes  is  accounted  for  by  multiplying  the  critical  velocity  ( Ucrk )  by  a 
correction  function  ( £//2 ).  It  is  noted  that  the  Soulsby- van  Rijn  (Soulsby 

1997)  formulas  are  very  similar  to  the  van  Rijn  (1984a, b;  2007a, b)  except 
for  the  definition  of  the  effective  velocity  and  the  recalibration  of  the  bed¬ 
load  formula  coefficients  in  van  Rijn  (2007a).  The  proposed  changes  for 
multiple-sized  sediments  should  be  verified  with  measurements  or 
numerical  simulations  for  non-uniformly-sized  sediment  transport. 

Watanabe 

The  equilibrium  total-load  sediment  transport  rate  is  determined  by 
Watanabe  (1987)  as 


Qv  =  [fsrs  +  fbO--rs)\psAwatU 


T  —  T 

b  max  c; 


pg 


(2-94) 
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where: 

qt„  =  potential  total-load  transport  rate  [kg/m/s] 
r  =  fraction  of  suspended  load  defined  by  Equation  (2-51)  [-] 
r6max  =  combined  wave-current  maximum  shear  stress  [Pa] 

Tcr  =  critical  shear  stress  of  incipient  motion  [Pa] 

AWat  =  empirical  coefficient  typically  ranging  from  0.1  to  2.0  [-]. 

The  critical  shear  stress  is  determined  from  Equations  (2-100)  and 
(2-101).  The  maximum  bed  shear  stress  ( r6max )  is  calculated  as  (Soulsby 

1997) 


hmax  =  VK+TwC0S(p)2+(Tu,sin(p)2  (2-95) 

where  (p  is  the  angle  between  the  waves  and  current;  rh  is  the  mean  shear 
stress  due  to  waves  and  currents;  and  tw  is  the  wave-related  bed  shear 
stress  which  is  calculated  here  using  Equations  (2-12)  and  (2-13). 

The  fraction  of  suspended  sediment  ( rs  )  is  estimated  using  the  van  Rijn 

(2007a, b)  transport  equations  described  above.  Besides  being  needed  in 
the  total-load  transport  equation  (Equation  2-50),  it  also  allows  the 
application  of  the  bed-  and  suspended-load  scaling  factors  in  a  way  similar 
to  all  other  transport  formula. 

The  Watanabe  (1987)  transport  formula  is  modified  for  multiple-sized 
sediments  as 


Qtk*  =  [fsrsk  +  f,  (!  -  rsk  )\psPuAvatU 


T 


b  max 


-  T 

?k  crk 


pg 


(2-96) 


where  the  subscript  k  indicates  the  sediment  size  class. 

Hiding  and  Exposure 

When  the  bed  material  is  composed  of  multiple  grain  sizes,  larger  grains 
have  a  greater  probability  of  being  exposed  to  the  flow  while  smaller  parti¬ 
cles  have  a  greater  probability  of  being  hidden  from  the  flow.  Figure  2-6 
shows  an  example  of  a  sediment  grain  (  d.)  being  hidden  by  dk . 
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Figure  2-6.  Schematic  of  the  exposure 
height  of  bed  sediment  grains. 


For  the  transport  formulas  described  above,  the  hiding  and  exposure 
mechanism  is  considered  by  correcting  the  critical  shear  stress  or  velocity 
using  a  hiding  and  exposure  correction  function  ( 4 ).  For  the  Lund-CIRP 

transport  formula,  an  alternate  approach  is  required  due  to  the  way  in 
which  the  Shields  number  and  grain  size  are  included  in  the  formulation; 
thus,  the  hiding  and  exposure  correction  function  is  directly  used  to 
multiply  the  transport  rate.  Two  methods  are  used  to  calculate  4  depen¬ 
ding  on  whether  the  sediment  transport  model  is  run  with  a  single  sediment 
size  or  with  multiple  sediment  sizes;  the  methods  are  described  in  the 
following  sections. 

Single-sized  sediment  transport 

In  some  applications,  the  coastal  bed  material  is  dominated  by  a  single 
sediment  size  with  patches  of  other  sediment  sizes  or  materials  (e.g.,  shell 
hash)  that  may  not  contribute  significantly  to  morphology  change  in  the 
areas  of  interest;  however,  they  may  modify  the  sediment  transport 
through  hiding  and  exposure.  For  example,  it  is  possible  for  the  bed 
material  to  consist  of  mostly  uniform  sand  with  patches  of  shell  fragments 
(bimodal  distribution)  in  some  regions.  The  shell  material  is  difficult  to 
model  numerically  because  it  is  usually  poorly  sorted,  and  its  hydraulic 
properties  are  unknown.  For  such  regions,  sediment  transport  models 
often  tend  to  over-estimate  erosion  since  the  impacts  of  hiding  effect  of  the 
coarser  shell  material  are  not  represented  (e.g.,  Cayocca  2001).  A  better 
and  more  physically  plausible  approach  is  to  use  the  local  bed  composition 
along  with  a  correction  to  account  for  the  hiding  and  exposure  effects  of 
the  uniform  sand  with  the  patches  of  coarser  shell  material.  For  single¬ 
sized  sediment  transport,  the  correction  function  for  hiding  and  exposure 
is  calculated  following  Parker  et  al.  (1982)  as 
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& 


(2-97) 


where  m  is  an  empirical  coefficient  between  0.5  to  1.0.  The 
aforementioned  sediment  transport  equations  are  implemented  by  using 
the  transport  grain  size  {dk )rather  than  the  bed  material  ( d50 ).  A  single 

and  constant  transport  size  ( dk )  is  used,  while  the  bed  material  ( d50 ) 
varies  spatially.  The  spatial  distribution  of  d50  can  be  obtained  from  field 

measurement  data  and,  for  simplicity,  is  assumed  constant  during  the 
model  simulation  time.  This  is  a  significant  assumption  and  may  not  be 
reasonable  for  some  applications.  However,  this  method  provides  a  simple 
conceptual  mechanism  for  considering  an  important  process  in  the 
proposed  single-sized  sediment  transport  model.  The  approach  has  been 
successfully  applied  to  Shinnecock  Inlet,  NY,  to  simulate  morphology 
change  at  a  coastal  inlet  (Sanchez  and  Wu  2011a).  A  more  accurate  and 
complex  approach  is  to  simulate  the  transport  and  sorting  of  multiple¬ 
sized  sediments. 


Multiple-sized  sediment  transport 

The  hiding  and  exposure  correction  for  each  sediment  size  class  is  based 
on  Wu  et  al.  (2000): 


& 


(2-98) 


where  m  is  an  empirical  coefficient  that  varies  for  each  transport  formula, 
approximately  equal  to  o. 6-1.0.  The  total  hiding  and  exposure  probabilities 
{Pek  and  Phk ,  respectively)  are  calculated  as 


N 

Phk=Y,Pu 

j= 1 


dk+dj 


N 

-  Pek=YJP 


^  ljdt  +  dj 


(2-99) 


where  N  is  the  number  of  grain-size  classes. 

Incipient  Motion 

In  the  case  of  the  Lund-CIRP  (Camenen  and  Larson  2005,  2007,  2008) 
and  Watanabe  (1987)  formulas,  the  incipient  motion  is  based  on  the 
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critical  Shields  parameter  and  estimated  using  the  formula  proposed  by 
Soulsby  (1997): 


0  s 

0  = - : - +  0.055  fl-  exp  (-0.02d) 

cr  1  +  1. 2d  L  v  ' 


(2-100) 


in  which  the  dimensionless  grain  size  (d)  is  defined  in  Equation  (2-66). 
The  critical  shear  stress  for  incipient  motion  is  given  by 


g(ps  ~  P)d 


@„ 


(2-101) 


The  critical  depth-averaged  velocity  for  currents  alone  ( Ucrc )  is  calculated 
using  the  formula  proposed  by  van  Rijn  (1984c): 


t/.  = 


0.19d50  log10 


8-5d°06log10 


Ah 


d 

Ah 


,  for  0.1  <  d50  <  0.5mm 


90  > 
\ 


90  ) 


(2-102) 


,  for  0.5  <  dn  <  2.0mm 


where  d50  and  d90  are  the  sediment  grain  size  in  meters  of  50th  and  90th 

percentiles,  respectively.  The  above  criteria  are  used  in  the  van  Rijn 
(2007a, b)  and  Soulsby-van  Rijn  (Soulsby  1997)  transport  formulas. 


The  critical  bottom  orbital  velocity  magnitude  for  waves  alone  is  calculated 
using  the  formulation  of  Komar  and  Miller  (1975): 


U 


crw 


0.24[(s-l)9f“d“sT“, 
0.95  [(s  - 1)9]"'  d™’T°" , 


for  0.1  <  dn  <  0.5  mm 
for  0.5  <  d50  <  2.0mm 


(2-103) 


where  Tp  is  the  peak  wave  period. 

Ripple  Dimensions 

The  bed  forms  calculated  by  CMS  are  the  wave-  and  current-related 
ripples.  The  ripple  height  (used  to  calculate  the  mixing  layer  thickness)  is 
estimated  as  the  maximum  of  the  current-  and  wave-related  ripple  heights 
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Hr  =  max(Hrc)Hr>J.  (2-104) 

The  current-related  ripple  height  and  length  are  calculated  as  (Soulsby 
1997) 


HrtC=Lrc/7  (2-105) 

Lrc  —  1000d5Q .  (2-106) 


The  wave-related  ripple  height  and  length  are  calculated  using  the 
expressions  proposed  by  van  Rijn  (1984b,  1989): 


H 


r,w 


0-22 Aw  for  ipw  <  10 

2.8x10  13(250-ipu,)5Aiu  for  10  <tyw  <  250 
0  for250<ipi(; 


(2-107) 


L 


r,w 


1.25 Aw  for  ipw  <  10 

1.4 x  10~6 (250  —  )2'5 Aw  for  10  <ipw  <250 

0  for  250  <ipw 


(2-108) 


where: 

lz  rr 

A  =  semi-orbital  excursion  =  [m/s] 

2  71 

u2 

ip  =  wave  mobility  parameter  = - - -  [-] 

(s-l)gd50 

d50  =  median  grain  size  [m] 

s  =  sediment  specific  gravity  [-] 
g  =  gravitational  constant  (9.81  m/s2) 

uw  =  bottom  orbital  velocity  [m/s]  (for  random  waves  uw  =  Jiurms ) 
T  =  wave  period  [s]  (for  random  waves  T  =  T  ). 

The  current-  and  wave-related  ripple  height  and  length  are  used  in 
calculating  the  bed  form  roughness  for  use  in  the  Lund-CIRP  transport 
formula. 
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Horizontal  Sediment  Mixing  Coefficient 

The  horizontal  sediment  mixing  coefficient  ( vs )  represents  the  combined 

effects  of  turbulent  diffusion  and  dispersion  due  to  non-uniform  vertical 
profiles.  In  CMS,  the  horizontal  sediment  mixing  coefficient  is  assumed  to 
be  proportional  to  the  total  eddy  viscosity  as 

ys=vt/(Js  (2-109) 

where  <js  is  the  Schmidt  number,  and  v,  is  the  total  eddy  viscosity.  There 

are  many  formulas  to  estimate  the  Schmidt  number;  however,  for 
simplicity,  it  is  assumed  to  be  constant  here.  The  default  value  for  the 
Schmidt  number  is  equal  to  1.0  but  may  be  modified  by  the  user. 

Boundary  Conditions 

At  closed  boundaries  (interface  between  wet  and  dry  cells),  the  sediment 
transport  rate  normal  to  the  boundary  is  set  to  zero.  The  inflow  boundary 
condition  requires  a  given  sediment  concentration  at  the  boundary. 
However,  for  most  coastal  applications,  the  actual  sediment  concentration 
is  not  available,  and  the  model  implements  the  equilibrium  concentration. 
The  equilibrium  concentration  is  the  concentration  that  is  reached  under 
steady  and  horizontally  uniform  conditions.  If  the  flow  is  directed  outward 
of  the  domain,  a  zero-gradient  boundary  condition  is  used  for  sediment 
concentration. 
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3  Numerical  Methods 

Overview 

CMS-Flow  has  both  implicit  and  explicit  solution  schemes.  The  explicit 
solver  is  designed  for  dynamic  problems  with  extensive  wetting  and  drying 
that  require  small  computational  time-steps,  while  the  implicit  solver  is 
intended  for  simulating  tidal-  and  wave-induced  circulation  at  tidal  inlets, 
navigation  channels,  and  adjacent  beaches.  A  detailed  description  of  the 
numerical  formulation  of  the  explicit  solver  of  CMS-Flow  can  be  found  in 
Buttolph  et  al.  (2006)  and  is  not  repeated  here.  The  sections  below 
specifically  refer  to  the  implicit  solver  of  CMS-Flow. 

The  implicit  solver  uses  the  SIMPLEC  (Semi-Implicit  Method  for  Pressure 
Linked  Equations  Consistent)  algorithm  (van  Doormal  and  Raithby  1984) 
on  a  non-staggered  grid  to  handle  the  coupling  of  water  level  and  velocity. 
Primary  variables  u-,  u-velocity,  and  water  level  are  stored  on  the  same  set 
of  grid  points,  and  fluxes  at  cell  faces  are  determined  using  a  Rhie  and 
Chow  (1983)  type  momentum  interpolation  method  (Wu  et  al.  2011).  The 
explicit  solver  uses  a  staggered  grid  with  velocities  at  the  cell  faces  and  the 
water  levels  and  water  depths  at  the  cell  centers  (Buttolph  et  al.  2006). 
CMS-Flow  also  calculates  salinity,  sediment  transport,  and  morphology 
change.  The  governing  equations  for  hydrodynamics  and  sediment  and 
salinity  transport  have  similar  forms  which  can  be  written  as  a  general 
transport  (advection-diffusion)  equation.  In  order  to  avoid  redundant 
derivations  of  discretized  equations,  the  discretization  of  the  general 
transport  equation  is  described  in  this  chapter,  and  the  same  discretization 
may  be  applied  to  all  of  the  transport  equations.  Then,  the  specific  solution 
procedures  for  hydrodynamics,  sediment  transport,  and  bed  change  are 
introduced. 

Computational  Grid 

The  explicit  time  scheme  in  CMS-Flow  supports  uniform  and  non- 
uniformly  spaced  Cartesian  grids  while  the  implicit  version  of  CMS-Flow 
supports  generic  Cartesian  grids  which  can  be  uniform,  non-uniform,  or 
telescoping.  The  telescoping  locally  refines  the  mesh  by  splitting  a  cell  into 
subcells.  The  only  requirement  imposed  by  the  numerical  methods  is  that 
the  cells  must  have  a  rectangular  shape.  Additional  requirements  are 
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imposed  by  the  user  interface  which  limits  the  variety  of  types  of  Cartesian 
grids  to  help  simplify  the  grid  generation  and  avoid  grid  quality  issues.  The 
following  requirements  are  applied: 

1.  Cells  can  only  be  subdivided  into  four  subcells. 

2.  Only  two  neighboring  cells  are  allowed  in  the  same  direction  (i.e.,  north, 
south,  east,  and  west). 

3.  Cells  may  have  a  maximum  of  six  neighbors. 

4.  Refinement  levels  must  be  spaced  by  at  least  one  cell  apart  (i.e.,  cells  that 
share  the  same  corner  must  be  one  refinement  level  apart). 

Examples  of  violations  of  the  above  four  requirements  are  shown  in 
sequential  order  in  Figure  3-1  from  a  to  d.  Limiting  the  number  of 
neighboring  cells  to  six  avoids  excessive  cell  refinement  and  limits  the 
band  width  of  the  coefficient  matrix  for  the  linearlized  system  of  equations 
solved  by  the  implicit  solution  scheme.  The  last  two  requirements  avoid 
having  excessive  cell  refinement  which  can  cause  numerical  instabilities. 
The  first  requirement  simplifies  the  grid  generation  process  but  may  be 
relaxed  in  future  versions.  The  last  requirement  is  enforced  for  grid  quality 
purposes  and  a  more  smoothly  varying  grid  resolution. 

Figure  3-1.  Examples  of  invalid  Cartesian  computational  grids. 


b. 
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The  CMS-Flow  grids  generated  in  the  SMS  interface  can  be  classified  as 
uniform,  non-uniform  Cartesian  grids,  regular  telescoping,  and  stretched 
telescoping  grids  (Figure  3-2). 


Figure  3-2.  Types  of  Cartesian  grids  supported  by  the  SMS  interface  and  CMS-Flow. 


One  important  aspect  of  incompressible  flow  models  is  the  location  of 
primary  variables:  velocity  and  pressure  (water  level).  On  a  staggered  grid, 
the  pressure  (water  level)  is  located  at  the  center  of  cells  and  the  u-  and  v- 
velocities  are  located  along  the  faces  of  cells  (Harlow  and  Welsh  1965; 
Patankar  1980).  On  a  non-staggered  grid,  all  of  the  primary  variables  are 
located  at  the  center  of  cells.  This  can  lead  to  the  spatial  oscillations  referred 
to  as  checkerboard  oscillations  and  model  instability  if  the  linear  interpola¬ 
tion  between  nodes  is  used  to  determine  the  fluxes  at  cell  faces.  A  staggered 
grid  can  more  easily  eliminate  these  oscillations  when  compared  to  a  non- 
staggered  grid;  however,  a  non-staggered  grid  involves  a  simpler  source 
code  and  can  minimize  the  number  of  coefficients  that  must  be  computed 
and  stored  during  a  simulation  because  many  of  the  terms  in  the  equations 
are  equal.  In  particular,  a  staggered  grid  is  more  complicated  when  handling 


ERDC/CHL  TR-14-2 


53 


the  interface  between  coarse  and  fine  cells  where  5-  or  6-face  control 
volumes  are  used.  Therefore,  a  non-staggered  (collocated)  grid  approach  is 
adopted  for  CMS-Flow  with  a  Rhie  and  Chow  (1983)  momentum 
interpolation  technique  used  to  eliminate  the  checkerboard  oscillations. 
Figure  3-3  shows  the  location  of  primary  variables  on  the  5-  and  7-point 
stencils  (computational  molecule). 


Figure  3-3.  Computational  stencils  and  control  volumes  for  two  types  of 
Cartesian  grids:  non-uniform  Cartesian  (left)  and  telescoping  grid  (right). 
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The  data  structure  for  a  grid  can  be  approached  in  three  ways:  1)  block- 
structured,  2)  hierarchical  tree,  and  3)  unstructured.  The  block-structured 
approach  divides  the  domain  into  multiple  blocks,  and  each  block  is 
treated  as  structured.  The  block  structured  approach  requires  a  special 
treatment  between  blocks  to  ensure  mass  and  momentum  balance  using 
this  approach.  The  hierarchical  tree  approach  is  memory  intensive  and 
requires  parent-child  relationships  and  a  tree  traverse  to  determine  mesh 
connectivity.  For  the  unstructured  approach,  all  cells  are  numbered  in  a  lD 
sequence,  and  tables  are  used  to  determine  the  connectivity  of  neighboring 
cells.  Among  these  three  approaches,  the  unstructured  approach  is  the 
most  simple  and  is  therefore  applied  in  CMS-Flow.  Computational  cells 
are  numbered  in  an  unstructured  manner  via  a  lD  index  array.  Inactive 
cells  (permanently  dry)  are  not  included  in  the  lD  index  array  to  save 
memory  and  computational  time.  All  active  computational  cells  are 
numbered  sequentially.  It  is  noted  that  using  an  unstructured  data 
approach  does  not  limit  the  model  from  using  matrix  solvers  designed  for 
structured  grids.  For  convenience  with  handling  boundary  conditions, 
each  boundary  cell  has  a  neighboring  ghost  cell  outside  of  the 
computational  domain. 
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General  Transport  Equation 

The  general  transport  equation  is  given  by 


d 

hip 

dt 

iP) 

+ 


d(hVjt/>) 

()X; 


Temporal  Term  Advection  Term 


d 

dxj 


Th 


dip 

dxj, 


+  § 

Source  Term 


Diffusion  Term 


(3-1) 


where: 

=  general  scalar 

h  =  total  water  depth  [m] 

P  =  correction  factor  [-] 

Vj  =  total  flux  velocity  [m/s] 

T  =  diffusion  coefficient  for  ip  [m2/s] 

S  =  source/sink  term  which  includes  all  remaining  terms. 

Note  that  in  the  case  of  the  continuity  and  momentum  equations,  (p  is  equal 
to  l  and  Vt ,  respectively.  In  the  case  of  sediment  and  salinity  transport,  tp  is 
equal  to  Ctk  and  Csal ,  respectively.  The  correction  factor  {ft)  is  equal  to  the 
total-load  correction  factor  and  equal  to  l  for  all  other  equations. 


Spatial  Discretization 

A  control -volume  technique  is  used  in  which  the  governing  equations  are 
integrated  over  a  control  volume  to  obtain  an  algebraic  equation  that  can 
be  solved  numerically.  Integration  of  Equation  (3-1)  over  a  control  volume 
(shaded  areas  in  Figure  3-3)  yields  the  following: 


r  d 

hip 

dA+f  9 

{  dt 

A 

[Pi 

A  UXj 

hV^-Thpp 

dx- 


dA=fsdA 


(3-2) 


d 

hp(pp 

dt 

Pp 

AAP  +^‘/z[(hI.V;.)^-r(h!.V^)]dL  =  SPAAP  (3-3) 


d 

hp(pp 

dt 

Pp  } 

Vf<Pf-Tf(VJ>)  ddf=SpAAp  (3-4) 
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where: 


A  f  A  A  -v 

ni={n1,n2) 

Vf=(niVi)f 

(f>P 

hp 

hf 

(H 

r/ 

A Ap 

A  If 
SP 


outward  unit  vector  normal  to  cell  face /[-] 

outward  cell  face  velocity  [m/s] 

value  of  (j)  at  the  cell  centroid 

total  water  depth  at  the  cell  centroid  [m] 

interpolated  total  water  depth  at  cell  face /  [m] 

outward  normal  gradient  of  /  at  cell  face/;  (V±^)f  =(njVi(/>)f 

interpolated  diffusion  coefficient  for  (/> 
area  of  cell  P  [m/s] 
length  of  cell  face /[m] 

source/sink  term. 


In  the  above  equations,  the  Gauss  Divergence  Theorem  has  been  used  to 
convert  the  area  integral  to  a  boundary  integral.  In  addition,  the  area 
integrals  have  been  evaluated  assuming  variables  vary  linearly  within  cells. 
The  cell  face  velocity  ( Vf )  is  calculated  using  a  momentum  interpolation 

method  similar  to  that  of  Rhie  and  Chow  (1983)  and  is  described  in 
Hych'odynamics  of  the  present  Chapter. 


Temporal  Discretization 

The  general  transport  equation  is  rewritten  as 


d 

hcj) 

dt 

[P) 

G 


(3-5) 


where  G  includes  all  the  remaining  terms.  For  stability  and  efficiency,  a 
fully  implicit  time-stepping  scheme  is  used  in  the  form 


1_ 

At 


(1  +  0.50) 


hn+1<f>n+1 

Pnl 


-  hnd)n 

(1  +  0)— 7/— +  0.50  ^  9 


Pn 


P 


n- 1 


Gn+1  (3-6) 


where  0  is  a  weighting  factor  between  o  and  1.  For  6=0,  the  scheme 
reduces  to  the  first-order  backward  Euler  scheme,  and  with  9  =  1,  the 
scheme  reduces  to  the  second-order  backward  scheme  (Ferziger  and  Peric 
1997).  The  superscripts  indicate  the  time-step  levels,  with  n+ 1  being  the 
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next  time-step,  n  being  the  current  time-step,  and  n- 1  being  the  previous 
time-step. 

Cell-face  Interpolation 

The  general  formula  for  estimating  the  cell-face  value  ( )  is  given  by 

4  =  A& +d-A)4+A(Wl +Ci-A)('i,V),  (3'7> 

where  fL  is  a  linear  interpolation  factor,  Vy  is  the  gradient  operator  in  the 
direction  parallel  to  face  f  and  y  is  the  distance  from  the  cell  center  to  the 
ghost  point  (O)  parallel  to  the  cell  face  ( / )  (Figure  3-4). 


Figure  3-4.  Schematic  showing  interface  interpolation:  (a)  coarse  to  fine  and  (b) 

fine  to  coarse. 


The  subscripts  ||  and  _L  indicate  the  components  parallel  and  normal  to 
the  cell  face.  By  definition,  ||=2|n1|-|-l|n2|  and  _L=l|n1|  +  2|n2|.Notethat 
for  neighboring  cells  without  any  refinement,  y  is  equal  to  zero,  and  the 

above  equation  is  consistent  with  non-refined  cell  faces.  The  linear 
interpolation  factor  is  defined  as 


f± 


x^f~x±,p 

X1,N  ~~  X±,P 


(3-8) 


where  x±  f  is  the  coordinate  of/perpendicular  to  the  face,  and  Axy  is  the 
cell  dimension  perpendicular  to  the  face/. 
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Cell-face  Gradient 

A  linearly  exact  second-order  approximation  for  the  normal  gradient  at 
cell  face /is  calculated  using  the  auxiliary  node  concept  of  Ferziger  and 
Peric  (1997): 


1  lW/"  N  W 


(3-9) 


where  the  subscripts  P  and  N  refer  to  two  neighboring  cells, 


<5  = 


X±,N  X±,P 


is  the  distance  between  cells  P  and  N,  normal  to  the  cell 


face  (Figure  3-4),  and  V,,  is  the  gradient  operator  in  the  direction  parallel  to 


face/.  Ham  et  al.  (2002)  compared  the  auxiliary  node  formulation  to  the 
fully-unstructured  discretization  proposed  by  Zwart  et  al.  (1998)  for  the 
viscous  terms  and  found  that  the  auxiliary  node  formulation  is  significantly 
more  stable. 


Cell-centered  Gradient 

The  cell-centered  gradient  operator  is  calculated  using  the  Gauss’ 
Divergence  Theorem  as 


J  V1/dA  =  J]hI4A//  (3-10) 

A  f 

The  equation  above  is  second  order  and  conservative  for  regular  and  non- 
uniform  grids. 

Reconstruction,  Monotonicity,  and  Slope  Limiters 

Variables  are  linearly  reconstructed  within  the  cells  as 

(f>  =  (f>p  +  r  -V  (fp  (3-11) 

where  /,  is  the  cell-average  value  specified  at  the  cell  centroid,  f  is  the 

distance  vector  from  the  cell  centroid  to  any  location  within  the  cell,  and 
V  (j)p  is  the  cell-centered  gradient.  The  reconstruction  is  second  order  and 

conservative  in  the  sense  that  <feAAp  =  /  <j>dA .  The  linear  reconstruction  is 
used  when  interpolating  cell-face  values  (Equation  3-7)  and  calculating 
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cell-face  gradients  (Equation  3-13).  If  the  reconstruction  satisfies  the  local 
maximum  principle 

min(^-^,0)<r  -V^  <max(^-^,0)  (3-12) 

then  no  new  extrema  are  created  within  the  cell,  and  the  solution  is 
monotonic.  Figure  3-5  shows  two  examples  of  linear  reconstruction  with 
and  without  slope  limiters  to  ensure  monotonicity. 


Figure  3-5.  Schematics  showing  examples  of  (a)  non-limited  and  (b)  limited  linear 

reconstructions. 


For  non-telescoping  grids,  the  slope  limiter  is  calculated  as 


<W)  = 


4  R 


(R+iy 

2  R 


(R2+ 1) 


van  Leer  (1979) 
van  Albadaetal. (1982) 


mm 


1, 


4  R 


R+l  R+l 


MUSCLfvan  Leer  1979) 


(3-13) 


where  R  is  the  ratio  between  two  consecutive  slopes 


R 


#+i -#)(*,■ 

xi+i~xi)(</{-<ti-i 


(3-14) 


Here,  the  second-order  van  Leer  (1979)  limiter  is  used  because  of  its 
smoothness.  Figure  3-6  shows  a  comparison  of  three  different  common 
limiters.  The  slope  limiter  is  applied  in  each  direction  separately. 
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Figure  3-6.  Comparison  of  three  different  slope  limiters. 


For  unstructured  grids  the  slope  limiters  described  above  are  difficult  to 
implement  because  of  the  difficulty  in  defining  forward  and  backward 
differences.  For  telescoping  grids  the  following  Limited  Central  Difference 
(LCD)  slope  limiting  procedure  of  Hubbard  (1999)  is  applied: 


max(4  -4,0) 
(r-V4P 
min(4— 4,0) 

(r-V4p 


1 


for  (r  •  V^)p  >max(4  0) 

for  (r  •  V  (ji)p  <  min(4  -  4  >  °) 
otherwise 


(3-15) 


where  rp=xf-xp.  In  the  procedure  outlined  by  Hubbard  (1999),  a  scalar 
limiter  is  then  calculated  as  <5  =  min^®^ .  Here,  a  directional  limiter  is 
applied  as  ®;  =  min  (®/ j ,  which  is  less  dissipative.  Finally,  the  cell- 
centered  gradient  is  limited  as 

vi4  =  ®iv;4  (3-16) 


_ _ . 

where  v;4  is  the  unlimited  gradient. 


ERDC/CHL  TR-14-2 


60 


Advection  Schemes 

Hybrid  Scheme 

The  hybrid  scheme  is  composed  of  a  first-order  upwind  scheme  and  a 
second-order  central  difference  scheme.  The  cell-face  advective  value  is 
given  by 


4  = 


(4  +</k)/ 2  for  \Pf\<2 
fa  for  |-f/|  >  2 


(3-17) 


where  the  subscripts  D  and  C  indicate  the  downstream  and  upstream 
values,  and  Pf  =  V f  |<5  /  Tf  is  the  Peclet  number  at  the  cell  face  in  which 


<5  = 


•Vl  ,n  x±,p 


.  In  the  hybrid  scheme,  when  the  Peclet  number  is  larger 


than  2,  the  first-order  upwind  value  is  used;  otherwise,  the  second-order 
central  difference  value  is  used. 


Exponential  Scheme 

The  exponential  scheme  interpolates  the  face  value  using  an  exact  solution 
to  the  lD  steady  advection-diffusion  equation: 


4~4  _exp(P//±)-l 
4  “4  exp(Pf )  —  1 


The  exponential  scheme  has  automatic  upwinding  and  is  stable  but  is 
usually  less  than  second  order. 

Hybrid  Linear/Parabolic  Scheme 

The  Hybrid  Linear/Parabolic  Approximation  (HLPA)  scheme  of  Zhu 
(1991)  may  be  written  as 


4 +(4-4)4  for  0  <4  <  1 

4  otherwise 


(3-19) 


where  the  subscripts  D ,  C ,  and  U  indicate  the  first  downstream  and  first 
and  second  upstream  cells,  respectively.  The  normalized  variable  (  4 )  is 
determined  based  on  the  formulation  of  Jasak  et  al.  (1999) 
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2  _  0c  0b  _  i _ 0p  0c 

2(V4§,C 


(3-20) 


where  8±c  —  x1D  —  x± >c .  The  HLPA  scheme  is  second  order.  Choi  et  al. 

(1995)  found  that  the  HLPA  scheme  has  similar  accuracy  to  the  third- 
order  SMARTER  (Sharp  and  Monotonic  Algorithm  for  Realistic  Transport 
Efficiently  Revised)  and  LPPA  (Linear  and  Piecewise-Parabolic 
Approximation)  schemes  but  is  simpler  and  more  efficient  (Shin  and  Choi 
1992;  Choi  et  al.  1995). 

Source/Sink  Term 

The  source/sink  term  is  linearized  as  (Patankar  1980) 

Sp=S£+S;4  (3-21) 


The  coefficient  Sp  is  required  to  be  non-positive  for  stability. 


Assembly  of  Algebraic  Equations 

Assembly  refers  to  the  process  of  combining  terms  to  create  a  linear 
system  of  algebraic  equations.  The  algebraic  equation  for  each  cell  is 
obtained  by  first  combining  or  assembling  all  of  the  terms.  Then,  the 
continuity  equation  is  multiplied  by  (/>"''  and  is  subtracted  from  the 

transport  equation.  The  resulting  discretized  equation  for  cell  P  is 

ap0?+1  =  EMT  +  b  {3-22) 

N 


where  the  subscript  N  refers  to  the  neighboring  cell  sharing  cell  faces;  ap 
and  aN  are  linear  coefficients  for  ^”+1  and  0”+1 .  The  last  term  ( b  )  contains 

all  the  remaining  terms.  Applying  a  similar  equation  for  all  of  the  internal 
cells  of  a  grid  results  in  a  system  of  algebraic  equations.  This  set  of 
equations  is  referred  to  as  the  discretized  governing  equations. 

Implicit  Relaxation 

The  right-hand  side  of  the  algebraic  system  of  equations  ( b )  generally 
contains  deferred  corrections  and  source  terms  which  must  be  computed 
iteratively.  This  iteration  loop  is  referred  to  as  the  outer  loop,  and  the  loop 
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within  the  linear  iterative  solver  (see  Section  Iterative  Solvers )  is  referred 
to  as  the  inner  loop.  Under-relaxation  is  often  applied  to  stabilize  the 
convergence  of  the  outer  iteration  loop.  Under-relaxation  maybe  applied 
as  (f>m+1  =  a^n+1  +(1  —  a^m  in  which  a *  is  an  under-relaxation 
parameter  (0  <a  <  1).  The  superscript  m  indicates  the  previous  iteration, 

and  superscript  m+ 1  indicates  the  new  relaxed  value.  A  better  approach  is 
used  here  to  introduce  the  under-relaxation  directly  in  the  algebraic 
system  of  equations  (Majumdar  1988;  Ferziger  and  Peric  1997)  as 

=  £a„$+1  +  6  +  ^^M,m  (3-23) 

S  N 

This  is  referred  to  implicit  under-relaxation  and  has  the  effect  of  making 
the  coefficient  matrix  more  diagonally  dominant,  thus  improving 
convergence. 

Iterative  Solvers 

The  selection  of  an  iterative  solver  is  one  of  the  key  issues  impacting  the 
overall  performance  of  the  model.  CMS  has  six  iterative  solvers  available, 
and  they  are  described  in  more  detail  below:  1)  GMRES,  2)  BiCGStab,  3) 

SIP,  4)  ICCG,  5)  Gauss-Seidel,  and  6)  Gauss-Seidel  with  Successive-Over- 
Relaxation.  The  default  iterative  solver  is  a  variation  of  the  GMRES 
(Generalized  Minimum  RESidual)  method  (Saad  1993).  The  original 
GMRES  method  (Saad  and  Schultz  1986)  utilizes  the  Arnoldi  process  to 
reduce  the  coefficient  matrix  to  the  Hessenburg  form  and  minimizes  the 
norm  of  the  residual  vector  over  a  Krylov  subspace  at  each  iterative  step. 

The  variation  of  the  GMRES  method  used  here  allows  changes  in 
preconditioning  at  every  iteration  step  (Saad  1993).  The  Incomplete  Lower 
Upper  Factorization  ILUT  (Saad  1994)  is  used  as  the  preconditioner  to 
speed  up  convergence.  The  GMRES  solver  is  applicable  to  symmetric  and 
non-symmetric  matrices  and  leads  to  the  smallest  residual  for  a  fixed 
number  of  iterations.  However,  the  memory  requirements  and  computa¬ 
tional  costs  become  increasingly  expensive  for  larger  systems. 

The  BiCGStab  (BiConjugate  Gradient  Stabilized)  iterative  solver  is  also  a 
Krylov  subspace  solver  and  is  applicable  to  symmetric  and  non-symmetric 
matrices  (Saad  1996).  BiCGStab  also  uses  ILUT  as  a  preconditioner  (Saad 
1994).  The  BiCGStab  method  can  be  viewed  as  a  combination  of  the 
standard  Biconjugate  Gradient  solver  where  each  iterative  step  is  followed 
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by  a  restarted  GMRES  iterative  step.  One  advantage  of  the  BiCGStab 
iterative  solver  is  that  the  memory  requirements  are  constant  for  each 
iteration,  and  there  are  less  computational  costs  when  compared  to  the 
GMRES  method  (for  less  than  four  iterations). 


The  SIP  (Strongly  Implicit  Procedure)  iterative  solver  uses  an  Incomplete 
Lower  Upper  decomposition,  which  approximates  the  exact  Lower  Upper 
decomposition  (Stone  1968).  The  method  is  specifically  designed  for 
algebraic  systems  of  equations  derived  from  partial  differential  equations. 
The  implementation  here  is  for  a  5-point  stencil  and,  therefore,  only 
applies  to  nontelescoping  grids. 

The  ICCG  (Incomplete  Cholesky  preconditioned  Conjugate  Gradient) 
iterative  solver  is  applicable  to  symmetric  matrices  such  as  the  pressure 
correction  equation  (Ferziger  and  Peric  2002).  The  implementation  here  is 
also  for  a  5-point  stencil  and,  therefore,  can  only  be  applied  to 
nontelescoping  grids. 

The  simplest  iterative  solvers  implemented  here  are  the  point-implicit 
Gauss-Seidel  solvers  with  or  without  Successive-Over-Relaxation.  The 
Succesive-Over-Relaxation  may  speed  up  convergence  but  can  also  lead  to 
model  divergence  (Patankar  1980).  Even  though  the  Gauss-Seidel  method 
requires  more  iterations  for  convergence,  the  overall  efficiency  may  be 
higher  than  the  GMRES  and  BiCGStab  because  each  iteration  is  computa¬ 
tionally  inexpensive,  and  the  code  is  efficiently  parallelized.  However,  based 
on  experience  and  testing,  the  GMRES  and  BiCGStab  methods  are  usually 
more  robust  and  perform  better  for  large  time-steps. 


Convergence  and  Time-Stepping 

During  the  iterative  solution  process,  error  is  calculated  and  used  to 
determine  if  a  solution  has  converged,  diverged,  or  stalled  at  an  error 
below  a  set  tolerance  threshold.  An  estimate  of  the  error  in  solving  the 
general  algebraic  equation  is  given  by 


rp 


1 

aP 


,  N 


\ 

-ap<t£+1  +b 


(3-24) 


The  /2-norm  of  the  residual  errors  is  given  by 
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(3-25) 


Since  this  value  depends  on  the  total  number  of  cells,  the  final  statistic 
(referred  to  as  the  residual)  that  is  used  for  estimating  the  model 
convergence  is  obtained  by  dividing  the  Z2-norm  by  the  square  root  of  the 
number  of  cells: 


Rm 


(3-26) 


where  Nc  is  the  number  of  computational  cells.  Simply  stated,  R'"  is 

referred  to  as  the  normalized  residual  error,  and  the  superscript  refers  to 
the  iteration  number.  Also,  Rm  is  calculated  for  each  variable  that  is  solved 
at  each  iteration  step  of  the  solution  process.  Each  equation  has  default 
maximum  tolerances  for  determining  if  the  solution  has  converged, 
diverged,  or  stalled.  The  maximum  number  of  iterations  that  is  imposed  is 
set  equal  to  M.  A  minimum  of  five  iterations  is  required  for  the 
hydrodynamic  equations,  and  a  minimum  of  M/2  iterations  is  required  for 
the  sediment  transport  equations.  Table  3-1  lists  the  default  criteria  to 
determine  whether  the  iterative  solution  procedure  has  converged, 
requires  a  reduced  time-step,  or  has  diverged. 


Table  3-1.  Default  criteria  to  determine  whether  the  iterative  solution  procedure  has 
converged,  diverged,  or  requires  a  reduced  time-step. 


Variable 

Converged 

Reduce  Time-Step 

Diverged 

Current  velocity,  m/s 

If  Rm<  1x10-7 

If  Rm>l. 0x10-3 

If  Rm>l. 0x10-2 

or  |  Ui  |  >10.0 

Pressu  re-correcti  on , 
m2/s2 

If  Rm<  1x10-8 

If  Rm>  1.0x10-4 

If  Rm  >1.0x10-3 

or  |  p  |  >50.0 

Total-load 
concentration, 
kg/ m3 

If  Rm  <1x10-8 

None 

If  Rm  >1.0x10-3 

or  Ctk<0 

Salinity,  ppt 

If  Rm  <1x10-6 

None 

If  Csal  <0 

For  the  implicit  model,  the  time-steps  for  the  hydrodynamics,  sediment 
and  salinity  transport  are  the  same  in  order  to  avoid  mass  conservation 
problems  and  for  simplicity.  If  any  of  the  time-step  reduction  criteria  are 
met,  then  the  time-step  is  reduced  by  half  and  a  minimum  number  of  two 
time-steps  are  calculated  at  the  newly  reduced  time-step.  If  the  last  time- 
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step  converged  properly,  then  the  time-step  is  increased.  The  maximum 
time-step  allowed  is  equal  to  the  user-specified  initial  time-step. 


Ramp  Function 

For  most  coastal  applications,  the  model  is  initialized  from  a  cold  start , 
which  means  that  the  water  level  and  current  velocities  are  initially  set  to 
zero.  When  the  model  is  initialized  with  anything  other  than  zeros,  this  is 
referred  to  as  a  hot  start.  The  ramp  period  allows  the  model  to  slowly 
transition  from  the  specified  initial  conditions  without  shocking  the 
system.  The  ramp  function  is  defined  here  as 


1  1 

f Ramp  =~Z~ 


2  2 


:7rmm 


,1 


L  Ramp 


(3-27) 


where  t  is  the  simulation  time,  and  tRamp  is  the  ramp  duration.  The  ramp 

function  provides  a  smooth  function  for  transitioning  from  the  specified 
initial  conditions  and  is  plotted  in  Figure  3-7. 


Figure  3-7.  Ramp  function  used  in  CMS. 


0  0.2  0.4  0.6  0.8  1  1.2 

t/t 

ramp 


The  ramp  function  is  applied  to  the  model  forcing  conditions,  including 
the  wave  forcing,  surface  wind,  sediment  concentration  capacity,  and 
significant  wave  height,  by  direct  multiplication  of  these  parameters  by  the 
ramp  function  at  each  time-step  during  the  ramp  period.  Boundary 
conditions  are  also  slowly  transitioned  from  the  initial  condition  by  direct 
multiplication  of  the  boundary  conditions  by  the  ramp  function  at  each 
time-step  during  the  ramp  period.  The  ramp  period  is  usually  on  the  order 
of  a  few  hours  to  a  few  days  and  should  not  be  confused  with  the  spin-up 
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period.  The  spin-up  period  is  the  time  it  takes  for  the  effects  of  the  initial 
condition  to  disappear  and  can  be  much  longer  than  the  ramp  period. 

Hydrodynamics 

Coupling  of  Velocity  and  Water  Level  -  SIMPLEC  Algorithm 

The  governing  equations  are  solved  in  a  segregated  manner  in  which  each 
governing  equation  is  solved  separately  in  a  sequential  manner  within  an 
iteration  loop  in  order  to  obtain  a  converged  solution.  Coupling  between  the 
velocity  (momentum)  and  water  level  (continuity)  is  achieved  with  the 
SIMPLEC  algorithm  (van  Doormal  and  Raithby  1984).  The  main  difficulty 
in  solving  the  momentum  equations  is  that  the  water  level  is  not  known  a 
priori  and  must  be  calculated  as  part  of  the  solution.  The  solution  algorithm 
procedure  is  described  below.  First,  an  initial  pressure  p  =  pgrj'  is 

estimated  based  on  the  previous  time-step  value.  Then,  the  momentum 
equations  are  solved  for  the  corresponding  velocity: 


dihVD 

o(hv;v;) 

\  J  ~ 

d 

vhw} 

h  dp*  ( 

dt 

dxj 

Pxj 

1  dx , 

v  J  ) 

p  dxi 

where  S*  includes  all  the  remaining  terms.  The  cell  face  velocities  are 
calculated  with  a  Rhie  and  Chow  (1983)  type  interpolation  method: 


-a,-. 


/  \ 

AAP 

(h„  * 

P 

— V±p 

ap  J 

/ 

iP 

(3-29) 


where  H  =  V  +  av  — — — V  p  ,  av  is  the  implicit  relaxation  coefficient 

ap  p 

(set  to  0.8  by  default),  and  (  )f  denotes  linear  interpolation  (see  Implicit 

Relaxation  of  the  present  Chapter).  This  method  avoids  the  checkerboard 
oscillations  associated  with  the  collocated  grid.  It  is  noted  that  this 
approach  is  slightly  different  from  that  originally  proposed  by  Rhie  and 
Chow  (1983)  and  used  by  several  authors  such  as  Lai  (2010).  This  current 
approach  was  found  to  be  significantly  more  stable. 

Next,  the  velocity  (V)  and  pressure  corrections  ( p' )  are  defined  such  that 
both  the  momentum  and  continuity  equations  are  satisfied: 
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V”+1  =  V*  +  V(  ,  pn+1  =  p*  +p'  (3-30  a,b) 


In  the  SIMPLEC  algorithm,  the  velocity  correction  is  related  to  the 
pressure  correction  by 


V[=  -GVfp' 


(3-1) 


a,. 


where  G 


/z  A Af 
P  0-p 


a 


N 


Using  dh/dt  =  (dp/  8t )  /  (pg)  and  substituting  Ij',+1  =V*  +  V/  in  the 
continuity  equation  yields  the  semi-discrete  water  level  correction 
equation 


(1  +  O.50)(p*  +  p’ )  -  (1  +  9)pn  +  0.5  9pn-1 
pgAt 


hG 


.(3-32) 


Note  that  at  convergence,  p  ->  o,  pn+1  =  p * ,  V/'+1  =  V* ,  and  the  above 

equation  reduces  to  the  continuity  equation.  Once  the  pressure  correction 
equation  is  solved,  the  cell-centered  water  levels  and  current  velocities  are 
corrected.  The  cell  face  velocities  are  also  corrected  as  VJ+1  =  V*f  +  V'f  ,  in 

which  the  velocity  correction  is  given  by  V'f  =  —Gf  (V  pr)  f . 


Summary  of  the  SIMPLEC  Algorithm: 

1.  Guess  the  water  level  and  pressure  field  ( p  ). 

2.  Solve  the  momentum  equations  (Equation  3-28)  to  obtain  V* . 

3.  Use  the  Rhie  and  Chow  (1983)  type  momentum  to  determine  the  velocities 
and  fluxes  at  cell  faces  (Equation  3-29). 

4.  Solve  the  pressure  equation  (Equation  3-34)  to  obtain  p' . 

5.  Use  the  correction  equations  to  adjust  the  velocities  and  water  levels. 

6.  Treat  the  corrected  water  level  and  pressure  field  as  a  new  guess  and 
repeat  this  procedure  from  Step  2  until  convergence  is  achieved. 
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Wetting  and  Drying 

During  numerical  simulations  of  the  surface  water  flows  with  sloped 
beaches,  sand  bars,  and  islands,  the  land-water  interface  changes  with 
time.  This  means  that  it  is  possible  for  nodes  at  the  land-water  interface  to 
be  wet  or  dry  throughout  a  given  simulation.  In  CMS,  a  threshold  water 
depth  is  used  to  judge  drying  and  wetting.  If  the  depth  at  the  cell  center  is 
larger  than  the  threshold  value  (i.e.,  a  small  value  such  as  0.02  m  for  field 
cases),  then  the  node  is  considered  to  be  wet.  Similarly,  if  the  depth  at  the 
cell  center  is  smaller  than  the  threshold  value,  then  the  node  is  considered 
to  be  dry.  For  the  implicit  solver,  all  of  the  wet  and  dry  cells  are  included  in 
the  matrix  solver.  Dry  cells  are  assigned  with  a  zero  velocity. 

Salinity  Transport 

Transport  Equation 

The  salinity  transport  equation  is  discretized  using  the  methods  described 
in  the  previous  section  entitled  General  Ti'ansport  Equation  and  are  not 
repeated  here. 

Laplace  Equation 

The  option  is  provided  to  estimate  the  initial  salinity  based  on  the  solution 
of  a  2DH  Laplace  equation  with  given  boundary  conditions  and  scattered 
salinity  values.  The  Laplace  equation  is  discretized  similarly  to  the 
transport  equation  as 


E(VxC^)A//=0  (3-33) 

/ 

where: 

(V1Ca,)/  =  outward  normal  gradient  of  (/)  at  cell  face  (/) 

A l  f  =  length  of  cell  face  ( / )  [m]. 

Sediment  Transport  and  Morphology  Change 

The  so-called  semi-coupled  sediment  transport  model  proposed  by  Wu 
(2004)  is  adopted  here,  in  which  the  sediment  calculations  are  decoupled 
from  the  hydrodynamic  calculations,  but  the  sediment  transport,  bed 
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change,  and  bed  material  sorting  equations  are  simultaneously  solved  in  a 
coupled  form  at  each  time-step. 

Transport  Equations 

The  sediment  transport  equations  are  discretized  using  the  methods 
described  in  the  previous  section  entitled  General  Transport  Equation 
and  are  not  repeated  here. 


Bed  Change  Equations 

The  fractional  bed  change  equation  is  discretized  as 


A z, 


•f 'morph 


bk 


Ps^-P’j 


afoUcr-Q, 


+  S, 


bk 


(3-34) 


where: 


f morph  =  morphologic  acceleration  factor 


Cn+1 
° bk 
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Ds  =  bed  slope  coefficient  [-] 

(qhk)f  =  magnitude  of  the  fraction  bed  load  at  the  cell  face/ calculated 

from  previous  iteration  [kg/m/s] 

V  zh )  =  bed  slope  calculated  at  the  cell  face. 


The  purpose  of  the  morphologic  acceleration  factor  ( fmorph )  is  to  speed  up 
the  bed  change  so  that  the  simulation  time  ( tsim )  represents  approximately 
the  change  that  would  occur  in  t  h  =  fmorph  tsim .  This  factor  should  be 

used  with  caution  and  only  for  idealized  cases  or  time  periods  which  are 
periodic  (mainly  tidal).  If  time-varying  winds  or  waves  are  important 
processes  for  driving  sediment  transport,  then  it  is  recommended  to  use 
reduced  or  representative  wind  and  wave  conditions.  Since  the  CMS  runs 
relatively  fast,  it  is  generally  recommended  not  to  use  the  morphologic 
acceleration  factor  when  validating  the  sediment  transport  model  using 
hindcast  measurements.  The  morphologic  acceleration  factor  is  useful, 
however,  when  simulating  idealized  cases  or  analyzing  project  alternatives. 
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Bed  Material  Sorting 

The  bed  material  sorting  equation  (Equation  2-53)  is  discretized  as 


pT 


AZbk±SMkZ_AZ2Pl_ 

on+l 


(3-35) 


where  Az2  =  Azb  -  <5'!+1  +  5"  is  the  change  in  the  top  elevation  of  the 
second  bed  layer,  and  p[n  =  p"k  for  Az2  >  0  and  p*kn  =  p2k  for  Az2  <  0 .  The 

mixing  layer  or  active  layer  thickness  calculation  is  slightly  modified  to 
avoid  excessively  small  layers  as 

5,  =  min[max(5lmin,2d50,A/2),5lmax]  (3-36) 


where  A  is  the  bed  form  height,  and  S^mi„  and  8l  max  are  the  user- 
specified  minimum  and  maximum  mixing  layer  thicknesses,  respectively. 

The  thickness  of  the  second  layer  is  calculated  as  82 +1  —  82  +  ISz2 .  The  bed 

material  gradation  in  the  second  layer  is  calculated  from  the  following 
discretized  form  of  Equation  2-54: 


n+1_8^pn2k+Az2pkn 
P2k  ~ 


8n+i 

u2 


(3-37) 


In  order  to  avoid  sediment  layers  from  becoming  extremely  thin  or  thick,  a 
layer  merging  and  splitting  algorithm  is  implemented  between  layers  2 
and  3.  Here,  the  subscript  2  corresponds  to  the  second  layer.  To  illustrate 
the  bed  layering  process,  Figure  3-8  shows  an  example  of  the  temporal 
evolution  of  seven  bed  layers  during  erosional  and  depositional  regimes. 

It  is  noted  that  in  order  to  maintain  a  constant  number  of  bed  layers,  the 
bottom  two  layers  are  merged  if  the  second  layer  is  split  or  a  new  layer  is 
created  at  the  bottom  using  the  bed  composition  and  thickness  of  the 
bottom  layer. 

Avalanching 

When  the  slope  of  a  non-cohesive  bed  ( <f)h )  is  larger  than  the  angle  of 
repose  {(/)R ),  the  bed  material  will  slide  (avalanche)  to  form  a  new  slope 
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Figure  3-8.  Schematic  showing  an  example  bed  layer  evolution.  Colors 
indicate  layer  number  and  not  bed  composition. 


approximately  equal  to  the  angle  of  repose.  The  process  of  avalanching  is 
simulated  by  enforcing  \<j)b\  <  <j>R  while  maintaining  mass  continuity 

between  adjacent  cells.  The  following  equation  for  bed  change  due  to 
avalanching  is  obtained  by  combining  the  equation  for  angle  of  repose  and 
the  continuity  equation  between  two  adjacent  cells: 

K,P  =  aa E  A  (tal1  ~  S§n  ^  tal1  ^  )H (1^  I  ~  ^  )  (3'38) 

where  SN  is  the  cell  center  distance  between  cells  P  and  N ,  AA  is  the  cell 
area,  aa  is  an  under-relaxation  factor  (approximately  0.25-0.5),  and  H(A') 

is  the  Heaviside  step-function  representing  the  activation  of  avalanching 
and  equal  to  1  for  X  >  o  and  o  for  X  <  o.  The  sign  function  ( sgn  X )  is  equal 

to  1  for  X  >  o  and  -1  for  X  <  o  and  accounts  for  the  fact  that  the  bed  slope 
may  have  a  negative  or  positive  sign.  Equation  (3-38)  is  applied  by  sweeping 
through  all  computational  cells  to  calculate  az;  and  then  modifying  the 
bathymetry  as  z'"-1  =  z'"  +  Az" .  Because  avalanching  between  two  cells  may 

induce  additional  avalanching  at  neighboring  cells,  the  above  sweeping 
process  is  repeated  until  avalanching  no  longer  occurs.  The  under¬ 
relaxation  factor  ( aa )  is  used  to  stabilize  the  avalanching  process  and  to 

avoid  overshooting  since  the  equation  is  derived  considering  only  two 
adjacent  cells  but  is  summed  over  all  (avalanching)  neighboring  cells. 
Equation  (3-38)  may  be  applied  to  any  grid  geometry  type  (i.e.,  triangles, 
rectangles,  etc.)  and  for  situations  in  which  neighboring  cells  are  joined  at 
corners  without  sharing  a  cell  face. 
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Hard  bottom 

The  sediment  transport  and  bed  change  equations  assume  a  loose  bottom 
in  which  the  bed  material  is  available  for  entrainment.  However,  hard 
bottoms  may  be  encountered  in  practical  engineering  applications  where 
bed  materials  are  non-erodible,  such  as  bare  rocks,  carbonate  reefs,  and 
concrete  coastal  structures.  Hard  bottom  cells  in  CMS  are  handled  by 
modifying  the  equilibrium  concentration  as  C't  =min (Ct„,Cf)  in  both  the 

sediment  transport  and  bed  change  equations.  The  bed  slope  term  in  the 
bed  change  equation  is  also  modified  so  that  only  deposition  (no  erosion) 
may  occur  at  hard-bottom  cells. 

Coupling  of  Sediment  Transport,  Bed  Change,  and  Sorting  Equations 

For  a  semi-coupled  sediment  transport  model,  the  sediment  calculations 
are  decoupled  from  the  hydrodynamics,  but  the  sediment  transport,  bed 
change,  and  bed  material  sorting  equations  are  coupled  at  the  time-step 
level  and  thus  solved  simultaneously.  A  modified  form  of  the  iteration 
procedure  of  Wu  (2004)  is  implemented  in  CMS.  The  equations  are 
obtained  by  substituting  Cff  =  p"+1C^n+1  into  the  bed  change  and  sorting 

equations  and  then  substituting  the  sorting  equation  into  the  bed  change 
equation. 

Summary  of  Sediment  Transport  Semi-Coupling  Procedure: 

1.  Calculate  bed  roughnesses  and  bed  shear  stresses. 

2.  Calculate  the  mixing  layer  thickness  S"+i . 

3.  Calculate  the  potential  sediment  concentration  capacity  C*f 1 . 

4.  Guess  the  new  bed  composition  as  Pik  =  Pik  ■ 

5.  Calculate  the  fractional  concentration  capacity  C"f  =  pf'Cff' . 

6.  Solve  transport  equations  for  each  sediment  size  class  for  C”f  . 

7.  Calculate  the  total  and  fractional  bed  changes. 

8.  Determine  the  bed  sorting  in  the  mixing  layer;  Repeat  Step  5  and  iterate 
until  convergence. 

9.  Update  the  bed  elevation  as  z"+1  =  znb  +  A zh . 

10.  Calculate  the  bed  gradation  in  the  bed  layers  below  the  mixing  layer. 

11.  Calculate  avalanching. 

12.  Correct  the  sediment  concentration  due  to  flow  depth  change. 
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When  using  the  explicit  time-stepping  scheme,  the  sediment  transport  and 
morphology  change  are  solved  at  a  time-step  which  may  be  equal  to  or 
multiples  of  the  hydrodynamic  time-step  for  efficiency.  However,  when 
using  the  implicit  time-stepping  scheme,  the  time-step  is  relatively  big  and 
on  the  order  of  to  min.  Because  the  time-step  is  so  big,  it  is  not  necessary 
to  use  different  time-steps  for  hydrodynamics  and  sediment  transport.  In 
addition,  using  the  same  time-step  for  hydrodynamics  and  sediment 
transport  provides  a  better  mass  balance.  For  these  reasons,  the  sediment 
transport  time-step  is  always  set  to  the  hydrodynamic  time-step  in  the 
implicit  time-stepping  scheme. 

Coupling  Procedure  of  CMS-Flow  and  CMS-Wave 

CMS-Flow  and  CMS-Wave  can  be  run  separately  or  coupled  together  using 
a  process  called  steering  (Figure  3-9).  The  variables  passed  from  CMS- 
Wave  to  CMS-Flow  are  the  significant  wave  height,  peak  wave  period, 
wave  direction,  wave  breaking  dissipation,  and  radiation  stress  gradients. 
CMS-Wave  uses  the  updated  bathymetry  (if  sediment  transport  is  turned 
on),  water  levels,  and  current  velocities  from  CMS-Flow. 


Figure  3-9.  Coupling  process  between  CMS-Flow  and  CMS-Wave. 


The  time  interval  at  which  the  CMS-Wave  model  is  run  is  called  the  steering 
interval.  Currently,  the  steering  interval  is  constant,  and,  therefore,  the 
input  spectra  must  be  at  constant  intervals  without  any  gaps.  A  schematic 
showing  the  steering  process  is  shown  in  Figure  3-10  in  which  the 
simulations  time  is  plotted  as  a  function  of  the  computational  time. 
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Figure  3-10.  Schematic  of  steering  process. 


QJ 

E 

i— 

c 

0 

4-> 

_03 

- CMS-Flow 

- CMS-Wave 

s 

1 

teering 

nterval 

/ 

E 

to 

/ 

Z 

w 

Computationa 

Time  r 

For  CMS  versions  prior  to  version  4.0,  the  steering  process  is  controlled  in 
the  SMS  interface  with  communication  files  because  the  CMS-Wave  and 
CMS-Flow  models  are  separate  executables.  For  CMS  V4.0  and  above,  both 
CMS-Flow  and  CMS-Wave  are  contained  within  a  single  executable 
(known  as  the  inline  executable ),  and  the  steering  process  is  controlled  by 
an  interval  steering  module.  Two  main  advantages  of  the  inline  executable 
and  inline  steering  module  are  as  follows:  1)  the  model  runs  faster  because 
there  is  no  need  to  use  communication  files  or  reinitialize  the  models 
(memory  allocation,  variable  initialization,  etc.);  and  2)  the  inline 
executable  makes  the  improvement  and  maintenance  of  the  steering 
module  easier  for  the  developers.  The  inline  steering  process  is 
summarized  as  follows: 

1.  CMS-Wave  is  run  for  the  first  two  time-steps,  and  the  wave  information  is 
passed  to  CMS-Flow  (Figure  3-10).  If  specified,  the  surface  roller  model  is 
run  on  the  wave  grid,  and  the  roller  contributions  to  the  radiation  stresses 
are  added  to  the  wave  radiation  stresses. 

2.  The  wave  height,  period,  dissipation,  radiation  stress  gradients,  and  wave 
unit  vectors  are  interpolated  spatially  from  the  CMS-Wave  grid  to  the 
CMS-Flow  grid. 

3.  CMS-Flow  is  run  until  the  next  steering  interval  and  wave  variables  are 
linearly  interpolated  throughout  time  during  the  specified  steering 
interval.  At  each  flow  time-step,  variables  such  as  wave  length  and  bottom 
orbital  velocities  are  updated  using  the  current  water  depths  and  current 
velocities. 
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4.  Water  levels,  current  velocities,  and  bed  elevations  are  estimated  for  the 
next  wave  time-step  and  are  interpolated  from  the  CMS-Flow  grid  to  the 
CMS-Wave  grid. 

5.  CMS-Wave  is  then  run  again  for  the  following  time-step. 

6.  Step  2-5  are  repeated  until  the  end  of  the  simulation. 

Spatial  Interpolation  and  Extrapolation 

CMS  allows  the  user  to  use  the  same  or  different  grids  for  CMS-Flow  and 
CMS-Wave.  If  the  same  grid  is  used,  then  no  spatial  interpolation  is 
carried  out.  If  different  grids  are  used,  then  spatial  interpolation  is 
necessary  in  order  to  transfer  information  from  one  model  grid  to  another 
model  grid.  The  interpolation  of  wave  variables  from  the  CMS-Wave  grid 
to  the  CMS-Flow  grid  is  done  using  a  combination  of  bilinear  and  linear 
triangular  interpolation.  Bilinear  interpolation  is  applied  at  non-jointed 
cells  (i.e.,  cells  with  four  neighbors)  and  triangular  interpolation  at  jointed 
cells  (cells  with  more  than  four  neighbors).  If  the  extents  of  the  CMS-Wave 
and  CMS-Flow  grids  are  different  (e.g.,  if  the  CMS-Flow  grid  is  smaller), 
then  the  extrapolation  of  variables  is  necessary  in  order  to  avoid  boundary 
problems  with  the  models.  The  spatial  extrapolations  for  different 
variables  are  given  by 


(3-39) 

(Ui’P  ^  wave  =  f  (Ui’N  )  flow 

(3-40) 

zw)L=(zw)H+/«(Az^)iw 

(3-41) 

(^s>p  ) flow  =  fex  K*  )wave 

(3-42) 

(t  y  ={t  y 

V  P’P  )  flow  V  P’^ /wave 

(3-43) 

( wp)n  =(wN)n 

V  p  /  flow  V  N  J  wave 

(3-44) 

where  ij  is  the  mean  water  surface  elevation,  Uj  is  depth-averaged  current 
velocities,  zb  is  the  bed  elevation,  Hs  is  the  significant  wave  height,  T  is 
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peak  wave  period,  and  w  is  the  wave  unit  vector.  The  superscripts  m  and  n 
indicate  the  CMS-Wave  and  CMS-Flow  time-steps,  respectively.  The 
subscripts  wave  and  flow  indicate  variables  on  CMS-Wave  and  CMS-Flow 
grids  respectively.  The  subscripts  P  and  N  indicate  variables  at  the 
extrapolated  and  nearest  neighbor  cells.  The  extrapolation  function  (  fext ) 

is  given  by 


1  1 

fext  =  2  +  2COS 


;rmin 

rN  y\ 

V 

ext 

(3-45) 


Here,  rN  is  the  distance  from  cell  P  to  N,  and  rext  is  the  extrapolation 

distance.  The  extrapolation  distance  is  assigned  to  each  computational 
grid  and  can  be  either  automatically  calculated  by  CMS  or  specified  by  the 
user.  The  mean  water  surface  elevation,  peak  wave  period,  and  wave  unit 
vectors  are  extrapolated  using  a  nearest  neighbor.  This  approach  is  more 
physically  accurate  than  extrapolating  only  to  a  certain  distance.  For 
example,  water  levels  are  controlled  mainly  by  tides  along  the  coast,  and 
the  spatial  variation  is  usually  much  smaller  than  the  tidal  range. 
Extrapolating  bed  elevations  from  a  boundary  can  lead  to  sharp  changes  in 
bathymetry  in  the  wave  model  and  instability  problems  in  both  the  wave 
and  flow  models.  A  better  approach  is  to  extrapolate  the  bed  change.  It  is 
noted  that  this  is  only  for  extrapolation.  For  internal  cells  the  actual  bed 
elevations  are  interpolated  from  the  flow  grid  to  the  wave  grid.  Finally, 
careful  attention  is  needed  in  determining  the  nearest  neighbor  so  that 
variables  are  not  extrapolated  over  inactive  portions  of  the  grid  (e.g., 
interpolating  values  in  a  bay  using  values  from  the  open  ocean). 

Temporal  Interpolation  and  Prediction 

Because  CMS-Wave  requires  the  water  surface  elevation  at  times  that  are 
ahead  of  the  hydrodynamic  model,  the  water  surface  elevation  and 
currents  must  be  predicted  for  the  CMS-Wave  time-step.  If  the  steering  is 
relatively  small  (<30  min),  then  the  values  from  the  last  time-step  may  be 
used  without  significant  error  as 


\  m 

ut) 

=  {uX 

(3-46) 

1  J flow 

\  1  /  flow 

(’7L 

={v)nflow 

(3-47) 
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where  superscripts  m  and  n  indicate  the  CMS-Wave  and  CMS-Flow  time- 
steps,  respectively.  The  subscript  flow  indicates  the  variables  on  the  CMS- 
Flow  grid.  For  many  coastal  engineering  applications,  it  is  desirable  and 
common  to  use  relatively  large  steering  intervals  of  2  to  3  hours.  For  large 
steering  intervals,  the  change  in  water  depth  has  the  largest  influence  on 
the  nearshore  wave  heights.  Therefore,  when  using  large  steering  intervals, 
it  is  desirable  to  make  better  predictions  of  water  levels  without  using  the 
previous  time-step  water  levels.  In  cases  where  the  relative  surface 
gradients  at  any  time  are  much  smaller  than  the  mean  tidal  elevation,  a 
better  approximation  of  water  level  may  be  obtained  by  decomposing  the 
water  level  into 


(3-49) 

where  rfm  is  the  mean  water  level,  and  ffv  is  a  variation  around  the  mean 

due  to  due  to  tidal,  wave,  and  wind  generated  surface  gradients.  The 
symbol  rjm  can  be  approximated  from  water  level  boundary  conditions  and 

is  generally  much  larger.  The  water  level  variation  ( rjv )  may  be  neglected 

or  approximated  using  the  last  flow  time-step  as 

feZ lw  ~  fe)^  =  fe^  "fe)i„,  (3-5°) 

For  most  coastal  inlet  applications,  the  above  equation  is  a  better 
representation  of  the  water  surface  elevation  and  is  used  as  the  default  in 
CMS.  After  spatially  interpolating  the  wave  variables  (significant  wave 
height,  wave  vector,  peak  wave  period,  wave  breaking  dissipation,  and 
wave  radiation  stresses)  onto  the  CMS-Flow  grid,  the  wave  variables  are 
linearly  interpolated  in  time.  The  wavelength,  bottom  orbital  velocities, 
and  mean  wave-current  bottom  friction  are  then  updated  for  current-wave 
interactions  at  each  flow  model  time-step. 
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4  Summary 

The  CMS  is  an  integrated  wave,  current,  sediment  transport,  and 
morphology  change  model  available  in  the  Surface-water  Modeling  System 
(SMS).  The  CMS  was  developed  under  the  Coastal  Inlets  Research 
Program  (CIRP),  an  Operations  and  Maintenance  Navigation  research 
program  at  the  Coastal  and  Hydraulics  Laboratory,  Engineer  Research  and 
Development  Center,  US  Army  Corps  of  Engineers.  The  CMS  consists  of 
two  main  components:  a  spectral  wave  transformation  model  named 
CMS-Wave  (Lin  et  al.  2008;  2011a, b)  and  a  hydrodynamic,  salinity, 
sediment,  and  morphology  change  model  named  CMS-Flow  (Wu  et  al. 
2010;  Sanchez  et  al.  2011a, b). 

Both  CMS-Wave  and  CMS-Flow  support  regular  and  non-uniform  Cartesian 
grids  while  the  CMS-Flow  supports  regular,  non-uniform,  telescoping,  and 
stretched  telescoping  Cartesian  grids.  The  CMS-Wave  and  CMS-Flow 
models  are  tightly  coupled  within  a  single  inline  code.  The  CMS-Wave  and 
CMS-Flow  grids  may  be  the  same  or  have  different  spatial  extents  and 
resolutions. 

CMS-Flow  uses  the  Finite  Volume  Method  and  has  both  fully  explicit  and 
fully  implicit  time-stepping  schemes.  Detailed  descriptions  of  the  CMS- 
Flow  mathematical  formulations  and  explicit  time-marching  schemes  are 
provided  in  Militello  et  al.  (2004)  and  Buttolph  et  al.  (2006)  and  are  not 
repeated  here.  Several  verification  and  validation  tests  were  presented  for 
hydrodynamics,  sediment  transport,  and  morphology  change  in  Sanchez 
et  al.  (2011a, b).  The  present  report  provides  an  updated  description  of  the 
mathematical  formulations  and  implicit  time  marching  schemes  which 
have  not  been  covered  in  previous  CMS-Flow  reports. 

The  hydrodynamic  model  solves  the  depth-integrated  and  wave-averaged 
hydrodynamic  equations  and  includes  physical  processes  such  as  advec- 
tion,  turbulent  mixing,  combined  wave-current  bottom  friction;  wave  mass 
flux;  wind,  atmospheric  pressure,  wave,  river,  and  tidal  forcing;  Coriolis 
force;  and  the  influence  of  coastal  structures  (Buttolph  et  al.  2006;  Wu 
et  al.  2010).  The  fully  implicit  hydrodynamic  model  uses  the  SIMPLEC 
(van  Doormaal  and  Raithby  1984)  algorithm  on  a  collocated  grid.  The 
inter-cell  velocities  are  calculated  using  a  Rhie  and  Chow  (1983)  type 
momentum  interpolation  method. 
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The  implicit  hydrodynamic  model  is  integrated  with  a  total-load 
nonequilibrium  and  multiple-sized  sediment  transport  model  and  includes 
processes  such  as  hiding  and  exposure,  bed  sorting  and  gradation,  bed 
slope  effects,  nonerodible  surfaces,  and  avalanching.  The  sediment 
transport  model  solves  a  depth -integrated  and  wave-averaged  total-load 
sediment  transport,  bed  sorting  and  change  equations  simultaneously 
using  an  iterative  method.  The  salinity  and  sediment  transport  equations 
are  solved  on  the  same  Cartesian  grid  as  the  hydrodynamics,  and  in  the 
case  of  the  implicit  time-stepping  scheme,  use  the  same  time-step  as  the 
hydrodynamics. 

The  latest  guidance,  documentation,  and  downloads  for  the  Coastal 
Modeling  System  are  available  at  http://cirp.usace.armv.mil/  and  from  the  CIRP 
wiki:  http://cirp.usace.armv.mil/wiki/Main  Page. 
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